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ABSTRACT 

Variational methods (VM) sensitivity analysis, which is the continuous alternative to the dis- 
crete sensitivity analysis, is employed to derive the costate (adjoint) equations, the transversality 
conditions, and the functional sensitivity derivatives. In the derivation of the sensitivity equa- 
tions, the variational methods use the generalized calculus of variations, in which the variable 
boundary is considered as the design function. The converged solution of the state equations 
together with the converged solution of the costate equations are integrated along the domain 
boundary to uniquely determine the functional sensitivity derivatives with respect to the design 
function. 

The determination of the sensitivity derivatives of the performance index or functional entails 
the coupled solutions of the state and costate equations. As the stable and converged numerical 
solution of the costate equations with their boundary conditions are a prioi unknown, numerical 
stability analysis is performed on both the state and costate equations. Thereafter, based on 
the amplification factors obtained by solving the generalized eigenvalue equations, the stability 
behavior of the costate equations is discussed and compared with the state (Euler) equations. 
The stability analysis of the costate equations suggests that the converged and stable solution 
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of the costate equation is possible only if the computational domain of the costate equations is 
transformed to take into account the reverse flow nature of the costate equations. 

The application of the variational methods to aerodynamic shape optimization problems is 
demonstrated for internal flow problems at supersonic Mach number range. The study shows, 
that while maintaining the accuracy of the functional sensitivity derivatives within the reasonable 
range for engineering prediction purposes, the variational methods show a substantial gain 
in computational efficiency, i.e., computer time and memory, when compared with the finite 
difference sensitivity analysis. 
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Chapter 1 
INTRODUCTION 


This chapter presents the global picture of design optimization and sensitivity 
analysis, briefly analyzes the historical and current state of optimization 
methodologies, outlines the motivations and objectives of the present research, 
and gives a short prelude to the remaining chapters . 

1.1 Overview of Aerodynamic Design Optimization and Sensitivity Analysis 

In the early times of flight, improvement of vehicle performance was mostly 
based first on intuition, empirically accumulated databases, and cut-and-try 
procedures [1,2]. Even recently, wind tunnel testing is being employed to perform 
optimization work to obtain airfoil performance criteria [3]. While this approach 
gave many valuable technical assistances, it was unable to furnish quick and 
reliable information to perform on-line design changes. 

In recent years, aerodynamic performance has been analyzed by a method of 
mathematical optimization. Optimizations can be performed either by those 
methods which need no gradient evaluations or by those which require gradient 
information. 

1.1.1 Nongradient Methods 

Variational methods (VM) were widely used to replace nonlinear partial 
differential flow equations and their boundary conditions with nonlinear algebraic 
equations and their corresponding boundary conditions so that approximate 
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solutions could easily be obtained. Bateman [4] is known for his formulation of 
the irrotational compressible inviscid flow using variational principles. Rasmussen 
and Heys [5] have extensively applied the application of variational methods to 
potential flows. 

Although the variational methods were used for the flow analysis in Refs. 4 
and 5, their application for design purposes is also well documented. With the 
Newtonian flow assumption and the linearized supersonic flow analysis, Miele [6 
-11] used the variational methods (VM) to obtain the optimality criteria. Then, 
employing this optimality criteria, he determined the geometry of a slender body 
of revolution having minimum pressure drag, optimized a two-dimensional wing 
for minimum pressure drag, designed optimal airfoils for supersonic speeds, and 
computed the optimal path for vehicles flying in different mediums. This approach 
has been used in solving many interesting and complex engineering problems 
with application to diversified fields, such as atmospheric and oceanographic 
studies [12,13] and planetary sciences [14]. To design nozzle shapes, Rao [15] 
combined variational approach and characteristic methods. Shmyglevskii [16] 
used VM and methods of characteristics to predict wave drag of high Mach 
number flow. Mikhlin [17] treated many mathematical-physics problems using this 
same method. In line with Rao [15], Thompson and Murthy [18] combined the 
characteristics methods and VM to design a three-dimensional rocket motor 
nozzle. In this class of optimization [6 - 18], the complete form of the Euler- 
Lagrangian equations and their boundary conditions are derived from the 
augmented Lagrangian and, thereafter, they are solved for the extremizing 
functions or curves until the measure of criteria is satisfied. 

From the combination of wide range of experiences in flow physics and wind 
tunnel techniques, the other category of optimization, i.e., the inverse design 
approach, has made immense contributions to design optimization. Lighthill [19], 
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in his pioneer work of optimization, showed how to numerically optimize a two- 
dimensinal airfoil for a prescribed pressure distribution using an inverse design 
approach. An ample number of researchers [20 - 24] have also applied this same 
method from lower subsonic to hypersonic flow regimes in two-and three- 
dimensional problems. Using a stream-line curvature method, Campbell [25] has 
also applied this approach to solve constrained optimization to design airfoils. 
The inherent problems with this approach, though, are the limitations in casting 
all relevant design problems in the form amenable to the inverse design 
optimization procedure and the requirement of high level of expertise to 
determine a priori the target objective functions. For this approach to work, the 
physics of the flow must be determined a priori in terms of the pressure or other 
quantities, and thereafter, the geometry that matches the above physical criterion 
must be sought. 

The other categories are the neural-network optimization approaches [26 - 
29]. While their applications are gaining momentum, they also have drawbacks 
due to the need for an extensive database, prior ideas about the optimal solution, 
intensive computations, and large computer memory. 

Note that all the design approaches mentioned previously do not require 
gradient computations. Therefore, they are efficient for moderate optimization 
problems. The limitations, though, are that those approaches are restricted to a 
certain class of problems with large databases and immense expertise and are 
confined to simplified geometries and flow field equations. 

1.1.2 Gradient-Based Methods (Numerical Design Optimization) 

With the advantage of modern hardware and software computer technologies, 
numerical design optimization and sensitivity analysis are currently solving 
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complete aircraft design in two-dimensional Navier-Stokes and three-dimensional 
Euler equations of the fluid flow. 

1.1 .2.1 Finite Difference Sensitivity Analysis 

The simplest, but the most expensive, sensitivity analysis technique used by 
gradient-based optimization methods is the finite difference approach. This 
method uses the one-sided or central-difference alternative to evaluate the 
sensitivities of performance functionals, and consequently, the computational 
time invested would increase with the increment of the number of design 
variables. This is due to the requirement to perturb each design variable by an 
appropriate step size and then compute the flow field variable for each new 
perturbed design variable with the chosen flow solver. This approach has an 
additional problem to determine the correct step size a priori so that the correct 
gradient is predicted within a given degree of accuracy. Despite its shortcomings, 
Huddelston and Mastin [30], and others, have applied this approach in their 
design procedure with Euler and Navier-Stokes equations as their flow field 
approximations. In the optimization package for general purposes optimizations, 
Vanderplaats [31] has also incorporated finite difference as an alternative to 
acquire the gradient information. 

1.1 .2.2 Discrete Sensitivity Analysis 

The other category of sensitivity analysis technique is the discrete analysis 
approach. The computation of the sensitivity equations is based on the Implicit 
Function Theorem. Due to the implicit dependence of the functional (objective 
function) and constraints on the flow field quantities, the determination of the 
sensitivity derivatives is related to obtaining the derivatives of the flow field vector 
with respect to the design variables. As the flow field equations are in most cases 
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solved in a computational domain, the functional dependence of the metric terms 
and the coordinate points with respect to the design variables are also required. 
This approach first calls upon the multiplication and assembly of various terms to 
a very large sparse linear algebraic equations, which depends on the number of 
design variables, and then solution of these sparse system of algebraic equations 
for the derivatives of the solution vector with respect to the design variables. 
Despite the large computational intensity and huge memory requirements of this 
approach, the versatility to incorporate many types of constraints, the need to 
perform multidisciplinary designs of moderate geometrical complexity, and the 
flexibility to incorporate it with any existing optimization algorithm make it 
attractive to perform design and shape optimizations. 

A wealth of literature can be found for this category. Hicks [32] and 
Vanderplaats [33,34] have used the discrete approach to design airfoils in 
transonic flow regimes. Pittman [35] has also used this procedure for supersonic 
flow conditions. Using the small perturbation equations in two dimensional flows, 
Elbana and Carlson [36] have also employed the technique. Recently Baysal and 
Eleshaky [37,38] and Eleshaky and Baysal [39] used this method for both internal 
and external flow problems. They also integrated the domain decomposition 
method in solving the sensitivity equations. Burgreen and Baysal [40,41] and 
Burgreen [42] further extended the methodology to the three-dimensional wing 
optimization and introduced an efficient way of parameterizing the curves and 
surfaces using the Bezier polynomials. Lacasse [43] applied the method to 
optimize multielement airfoils for two-dimensional, thin-layer, Navier-Stokes 
laminar equations. With a variant of approximation to the fluid flow, Taylor et al. 
[44,49], Newman et al. [45], Korivi et al. [46,48], and Hou et al. [47] introduced an 
incremental iterative technique to obtain the gradient information. In doing so, 
they have applied this new approach not only to the two-dimensional Euler and 
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thin-layer Navier-Stokes turbulent equations for internal and external flows but 
also to the three-dimensional Euler equations in supersonic flow regimes. 

1.1 .2.3 Variational Sensitivity Analysis 

The new emerging sensitivity analysis technique for gradient-based 
optimization methodology within the optimization community is the continuous 
sensitivity (variational sensitivity) analysis which fully exploits the variational 
methods. From the modified functional, this approach derives a set of partial 
differential equations (PDEs), i.e. the costate equations with their boundary 
conditions and the sensitivity equations. In computing the sensitivity derivatives 
with respect to the control points or design variables, this approach makes use of 
the converged solution of the state and costate equations. 

In recent years, variational sensitivity analysis has significantly contributed to 
the progress of aerodynamic design optimization. Lions [50], Pironneau [51], and 
Glowinski and Pironneau [52] showed the usefulness of the variational approach 
in fluid mechanical problems by illustrating how to compute the minimum drag 
profile in two-dimensional viscous and laminar flows. Chen and Seinfeld [53) 
developed a methodology to compute the performance sensitivity derivatives 
using optimal control theory. Koda et al. [54] used this procedure to solve 
atmospheric diffusion problems. Koda [55 - 57] further developed this approach 
and outlined a numerical algorithm for the computation of functional derivatives. 
This approach is well suited to solving the optimum design problems in fluid 
mechanics. Meric [58,59] treated optimal control problems governed by parabolic 
and elliptic partial differential equations and solved them numerically using 
variational methods. In their effort to compare the gradients obtained by "implicit" 
and "variational" approaches, Shubin and Frank [60] implemented VM to optimize 
the shape of a nozzle of a variable cross - sectional area for steady 
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one-dimensional Euler equations. Hou and Sheen [61] used a class of VM to 
derive second-order shape sensitivity equations of heat conduction problems. 
Jameson [62] regarded the boundary of the flow domain as a control parameter 
and then designed airfoils using the potential as well as the two- and three- 
dimensional compressible inviscid flows. Cabuk and Modi [63] implemented a 
perturbation method to compute the optimum profile of a diffuser for a maximum 
static pressure in a two-dimensional steady viscous incompressible flow. Ta'asan 
et al. [64] have successfully implemented variational methods and optimized an 
airfoil in the potential flow field. Quite recently, Ibrahim and Baysal [65] 
demonstrated the versatility of the variational methods to solve aerodynamical 
design problems for internal flows in different Mach number regimes including 
shock flows. Following the same approach as Jamson [62], Reuter and Jameson 
[66] optimized airfoils in potential flows. Also following the same solution method 
of Ta'asan et al. [64], Kurivila et al. [67] used the potential equations as their state 
equations and optimized the NACA 0012 airfoil for a given pressure distribution, 
lollo and Salas [68] used variational methods to solve two-dimensional internal 
flow optimization problem with embedded shock to match a pressure distribution. 
In addition to the general application of VM, many researchers have also 
proposed a numerical algorithm to accelerate convergence and improve 
efficiency for optimal design problems [69,70]. In this class of optimization, the 
functional sensitivity derivatives are directly coupled to the solution of a set of 
linear partial differential equations, i.e., the costate equations and their boundary 
or transversality conditions that result from the variation of the augmented 
Lagrangian function. The success of any optimization by this approach is, 
therefore, destined to a stable and converged solution of the costate equations. 
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1.2 Stability Analysis of Euler and Costate Equations 

The most popular schemes to advance any PDE, such as the Euler equations 
to steady-state solutions, are the implicitly factored time integration schemes 
(ADI). Unfortunately, approximate factorizations introduce errors that propagate 
throughout the computational domain. As a result of this, the stability limit is 
drastically reduced and the convergence rate deteriorates. To propose the range 
of Courant-Friedrichs-Lewy (CFL) numbers for which the allowable maximum 
eigenvalues are predicted, a stability analysis of the Euler and costate equations 
for the optimization purposes is conducted. 

Jesperson and Pulliam [71] studied the stability characteristics of the Euler 
equations for different flux-splitting methods. Anderson and Thomas [72] further 
conducted stability analysis on the complete three-dimensional Euler equations. 
Demuren and Ibraheem [73] have also pursued an extensive and complete 
stability analysis of one-, two-, three-dimensional Eluer and two- and three- 
dimensional Navier-Stokes equations. The common conclusion of these 
researchers [71 - 73] is that the stability solutions of state equations are impaired 
becauase of factorization and are dependent on the types of splitting and flux 
approximations. Also, the stability deteriorates as the dimensionality of the fluid 
equations increases. The stability problem associated with the numerical 
integrations of the costate equations was reported in Ibrahim and Baysal [65]. 
Although no stability analysis was performed for the costate equations in Ref. 65, 
its convergence and stability were assured by taking into account the reverse 
flow nature of the costate equations. 
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1.3 Motivations 

The advantage and disadvantage of the various approaches to perform 
design optimizations and sensitivity analysis are briefly discussed in Sec. 1.1. 
The main thrust of the exposition is intended to understand the complexities 
involved in each approach and to take advantage of the best of each approach 
as complementary to the others to realize a certain level of optimization goals. 

The major problems associated with any gradient-based numerical 
optimization are the penalties incurred because of the computational memory 
and time in obtaining the sensitivity gradients. These bottle-necks can be 
attributed to the repetitive need to solve the flow analysis for a change of any 
design variable and the huge memory allocation needed to store the derivatives 
of the objective function and constraints with respect to the design variables. This 
problem becomes more acute when one wants to design complex geometry 
using full Navier-Stokes equations at many design points in a multidisciplinary 
mode. While enjoying wide popularity and applicability, discrete sensitivity 
analysis has the previously mentioned potential limitations. Although no 
methodology exists which would overcome the aforementioned pitfalls of 
numerical optimization problems with the current state of computing facilities, 
sensitivity analysis by the variation methods (VM) is proposed in this study to 
partially alleviate the problems associated with huge memory allocation for 
moderate two-dimensional Navier-stokes and three-dimensional Euler equations. 
In this study, the motivation can be streamlined into three sub groups: efficiency 
in memory, efficiency in time, and generality in application. 

1.3.1 Computational Efficiency 

Independent of the approach one adopts, the complete optimization of 
gradient-based approach requires a repetitive solution of the state equations. 



10 


This being common to all, one is also concerned with the computational intensity 
and memory efficiency in the pursuits of the optimal solution. Here, efficiency is 
measured by the number of mathematical operational counts (time) associated 
with the computation and the computer memory saved in using the new proposed 
approach. The sensitivity analysis by the variational methods (VM) proposed in 
this study overcomes some of the critical issues by having some flexibility. First, 
the sensitivity analysis by VM involves the solution of the costate equations with 
their boundary conditions for already converged solution of the state equations. 
Also, because there are no metric and grid sensitivities as a result of a unique 
feature of the sensitivity analysis by variational methods, a substantial saving in 
computer memory is achieved. This is detrimental in large two-dimensional 
Navier-stokes and three-dimensional Euler equations. 

1.3.2 Generality of the Variational Approach 

First, since the costate equations are once and for all derived from the 
continuous PDE of the state equation, any robust solution method can be 
adopted to furnish the converged solution so that the costate equations can be 
solved until convergence is attained. This means that one does not necessarily 
have to solve the original state equation from which the costate equations are 
derived. Secondly, any other convenient discretization methods different from the 
type of discretization one uses for the state equations can be selected for the 
costate equations. The requirement that the costate equations be discretized 
exactly the same way as the state equations is shown not to be necessary, at 
least for quasi one-dimensional Euler equations [74], Thirdly, any time integration 
method different from the time integration method used for the state equation can 
be selected to advance the costate equations to steady state. The fourth point to 
mention is the design variables. In the approach proposed, note that the shape of 
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the domain is considered as the design parameter, and its contribution to the 
functional sensitivity derivatives is directly incorporated as shown in Chap. 3 and 
Appendices B and C. 


1.4 Objectives 

From the short motivations mentioned in Sec. 1.3, one would set the following 

objectives for the this study: 

(1) Develop a design optimization methodology based on sensitivity computed 
by the variational methods (VM), which is computationally efficient and 
general for aerodynamic optimization applications. 

(2) Investigate the causes of slow convergence and establish stability criteria 
for the costate equations by performing a complete stability analysis of the 
costate and Euler equations. 

(3) Demonstrate the concept on quasi one-dimensional Euler flow problems. 

(4) Extend the methodology to two-dimensional Euler equations of internal 
flows. 


1.5 Prelude to Chapters 2-7 

Chapter 2 presents the general two-dimensional Euler equations of the fluid 
flow in conservative forms as used in the general purpose CFL3D computer code 
[75] and the quasi one-dimensional Euler equations in conservative and 
nonconservative forms [76]. Chapter 3 gives a detailed procedure for 
aerodynamic sensitivity analysis of a two-dimensional nozzle problem by the 
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variational methods including the pertinent formulations and, finally, the 
derivation of the costate equations, their boundary conditions, and the sensitivity 
equations. Chapter 4 addresses the numerical approximations, discretizations, 
and time integrations of the costate equations. As the numerical solution of the 
costate equations are the prerequisite to obtain the gradient information, a 
complete numerical stability analysis is presented, and the results are discussed 
in Chap. 5 for the one- and two-dimensional equations of the costate and Euler 
equations. Chapter 6 demonstrates the applications of variational methods for 
sensitivity analysis by applying it to the quasi one-dimensional Euler equations 
for conservative and nonconservative flow field quantities. Chapter 7 discusses 
the optimization results for two-dimensional internal flow problem of the proposed 
approach for aerodynamic sensitivity analysis, and Chap. 8 finalizes the 
proposed theme (variational methods) and gives some recommendations for the 
future work in this particular area. 
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Chapter 2 

FLOW FIELD ANALYSIS 


2.1 Rationale for Two-Dimensional Euler Equations 

One of the integral parts of an optimization procedure is the correct 
approximation of the state equations. Depending on the simplicity, efficiency, and 
accuracy of the optimization procedure, any level of approximation to the flow 
physics can be adopted. To demonstrate the proposed optimization 
methodology, i.e., the variational methods (VM), we have chosen the Euler 
equations in conservative form. The choice of the state equations in conservative 
form renders it simpler to derive the costate equations along with their auxiliary 
conditions and sensitivity coefficients. 


2.2 Governing Equations of Two-Dimensional Euler Equations 

The Two-Dimensional, unsteady, Euler equations can be written in an integral 
form as [75,77] 



where G represents the vector of fluxes that is represented by 
G = E°i+F°j 


(2.2) 
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Here, E, and F are the cartesian flux components and i , and j are the unit 
vectors in the x and y directions, respectively. In Eq. (2.1), N is the unit normal 
pointing outward at the boundary and is defined as 

N = \n x ,n y ] r (2.3) 


The conservative dimensionalized solution vector and fluxes of the two- 
dimensional Euler equations in the Cartesian coordinate systems are given as 

Q = [p,pu,pv,pe,] T (2.4) 

E = [pu,pu 2 + p,puv,(pe, + p)u] T (2.5) 

F = [pv,pvu,pv 2 + p,(pe, + p)vf (2.6) 


where p, u, v, e t , p, E, and Fare the density, velocity in the x direction, 
velocity in the y direction, total internal energy, pressure, flux in the x direction, 
and, flux in the y direction, respectively. 

By an assumption of a smooth continuity on the integrand and application of 
the divergence theorem, Eq. (2.1) can be transformed to its mathematically 
equivalent differential form as 


+ ^ = 5 

dt dx dy 


(2.7) 


Given the reference length L, the free-stream density p„, and the free-stream 
speed of sound a„ the dimensional solution vector Q, and the other relevant 
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dimensional quantities in Eq. (2.7) are nondimensionalized with the help of the 
following definitions: 


x 




(2.8a) 



(2.8b) 


P = 


P_ 

P- 



(2.8c) 


where the tilda denotes a dimensionalized quantity. To perform time integration in 
uniform meshes, Eq. (2.7) is further transformed to a boundary-conforming space 
by use of stationary generalized curvilinear coordinate transformation as follows: 

£ = £(* ,y) Tl=V(x,y) x =t (2.9) 


With the help of Eqs. (2.8) and (2.9), Eq. (2.7) is transformed to 


dQ dE dF 
dt + dt; + dr\ 


= 0 


( 2 . 10 ) 


where 




Q = y = r'[p,pu,pv,pe$ 


( 2 . 11 ) 


( 2 . 12 ) 
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v 

(2.13a) 

V = T] x u+ T] y V 

(2.13b) 

E = 7 -1 [pU,puU + P£ x ,pvU + PZ y ,(P + e, )£/f 

(2.14a) 

F = J~ l \pV,puV + Pr\ x ,pvV + PT) y ,(P + e,)V] T 

(2.14b) 

In Eqs. (2.12) - (2.14), Q is the vector of conservative variables, E and F are the 
fluxes, U and V are the contravariant velocities in the new coordinate systems, 
and J is the Jacobian of transformation. The pressure P is related to the field 
variables through 

P = (r-\)p{e,-0.5(u 2 +v 2 )} 

(2.15) 

The various metric terms in Eqs. (2.9) - (2.14) are computed as 


= Jy n Vx = -Jy t = - Jx v % = Jx i 

(2.16) 


2.3 Solution Algorithm 

The basic governing equations and their formulations in the finite volume 
sense are presented in this section. Also, the various steps in approximating and 
advancing the solution to the steady state are briefly discussed. 

2.3.1 Finite Volume Formulation 

The guarantee to satisfy the mass, momentum, and energy across the cell 
faces, the ease with which to deal with complex geometries and discontinuities, 
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and the choice to work either on the physical or computational domain make the 
finite volume approach preferable to the finite difference approach. In this two- 
dimensional approach, the conserved field variables are cell-area averaged and 
computed at the cell center while the fluxes are computed at the cell faces with 
these cell-centered quantities. If one integrates Eq. (2.6) over the control volume 
bounded by lines of constant £ and 77 , then the resulting semidiscrete equations 
become: 


dQ, 

dt 


( + E ‘+y 2 .j E ‘-y 2 .i + F i.j*y 2 F ‘.j-y 2 ~ ® 


(2.17) 


where / and j are the grid points. We have also taken 


= 1 A77 = 1 


(2.18) 


and the cell-averaged Q t / is computed as 

tw} 


(2.19) 


2.3.2 Inviscid Upwind Spatial Differencing 

A Monotone Upstream Centered Scheme for Conservation Law (MUSCL) is 
generally adopted in the CFL3D computer code [75]. For instance, the derivative 
in the £ direction is 

(— ) = E —E 

K dt ) * i+ y* j 


( 2 . 20 ) 
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where the cell-face fluxes are constructed from the cell-center solution variables, 
which are located on the left and right sides of the cell faces. Mathematically, the 
fluxes can be expressed as 




( 2 . 21 ) 


where q+ are the nonconservative variables constructed from the upwind biased 
interpolations given by 




(2.22a) 


Ki = - ll[ (1 " * )A + + (1 + * )A -ll 

<+ 2’^ ^ ^ J i+1 J 


(2.22b) 


with k assuming three different values depending on the order of approximation. 
For instance, k = 1 for central difference, = 1/3 for the third-order upwind-biased, 
and = -1 for the second-order fully upwind. For flows with large flow field 
discontinuities, such as shocks, flux limiting of Chakravarthy-Osher [77] are used 
to maintain monotonocity. This can be achieved by 


A + = max[0,min(A + ^«A_,^A_jgrtA + )]^«A + 


(2.23a) 


A_ = max[0,min(A_5gnA +> /3A + ^g/jA_)]^«A. 


(2.23b) 


(3 ~£) 
(1 ~k) 


(2.24) 
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and 

q = [p,u,v,P] T (2.25) 

Also, a number of other upwind-biased interpolations are available in the 
literature, such as the Spekreijse, ENO, Venkat, and Superbee [78]. 

2.3.3 Van Leer Flux-Vector Splitting 

The most popular flux-vector-splitting schemes in Computational Fluid 
Dynamics (CFD) are the Beam-Warming [79] and Van Leer [80] flux-vector 
splittings. The Beam-Warming fluxes are constructed from the approximate 
fluxes. The Van Leer fluxes, on the other hand, are computed from the exact 
fluxes that fulfill definite criteria to maintain continuous differentiabilty at sonic and 
stagnation points. Because of this advantage over the Beam-Warming flux- 
vector-splitting method, the Van Leer method is adopted in this analysis. In 
practical computation, the Van Leer fluxes and flux Jacobians are split into 
positive and negative contributions by use of the local Mach number computed 
normal to the cell face. The local Mach number, in the £ direction for instance, is 
computed as 

u -m <2 - 26) 

where a is the local speed of sound and U is the contravariant velocity. For a 
supersonic flow where M^| > l, the positive and negative fluxes are expressed as 
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E + = E, E~ - 0 for M { > 1 (2.27a) 

E~ = E E + = 0 for < -1 (2.27b) 

For subsonic flow where < 1, the flux in the normal direction of the cell for 
constant t, is given by 

E = E + (Q-)+E~(Q + ) (2.28) 


The split fluxes in the 77 direction can be likewise obtained by replacing £ with 

Tl- 

From Eqs. (2.24) - (2.28), the split Jacobians in the £ and 77 directions, 
respectively, are obtained from the Van Leer fluxes, which are given in Appendix 
E, as 


dE\QT) 


dE~(Q + ) 


A-A + + A' 


(2.29) 


i, + _dF + {Q-) 
B 



B = B + + B~ 


(2.30) 


2.3.4 Time Integration 

Numerical integration in time can be accomplished by either explicit or implicit 
schemes. The choice of the type of time integration depends mostly on 
convergence, stability, efficiency with computer memory, and the flow physics. 
Implicit schemes are robust and stable for any type of flow and are not restrictive 
in the range of Courant-Friedrechs-Lewy (CFL) numbers. However, the explicit 
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schemes are inexpensive but restrictive with the allowable CFL numbers. To 
increase the range of the CFL number and assure stability, either dissipations are 
added or a sort of residual smoothing is performed [81]. In the case where the 
flow is unsteady, time-accurate explicit schemes are preferable even though 
implicit schemes with large time steps are often used as well. Note that, even for 
steady flow cases, the solution is marched in pseudo-time where the time is used 
as a relaxation factor. 

The application of backward Euler linearization in time whereby the higher- 
order terms are neglected, the flux E in the £ direction can be approximated as 


E H+l (Q) = E n (Q) + 



(2.31) 


Equation (2.31) can be simplified further as 


SE-(Q) bq 


E"‘(Q) = E'(Q) + —4^ ~Ar 
dQ dx 


(2.32) 


Now, the partial derivative of Q with respect to time in Eq. (2.32) is approximated 
by 


dQ e -<2" A<2 


dx Ax A-r 


Eq. (2.31) is simplified to 


, , _ . _ dE n (Q) _ 

E ( Q ) = E n (<2) + — jc 1 ■ AQ n 

dQ 


(2.33) 


(2.34) 
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where n is the time level. If one does the same type of linearization and 
approximations for the fluxes in the 77 direction and puts the approximate 
equations into Eq. (2.7), one obtains the following equations in the delta form: 

{/ + Ar[<5 { + A _ + S'A + + S;B + + <5;/r]}A{2 n = -At - *"((?) (2.35) 

where 

£"(g) = [S^t + 8-t + S;f- + <5;F + ] n (2.36) 


is the residual and 


A Q n = Q n+l -Q n (2.37) 

is the increment in the solution vector at each time level. 

For an initially guessed solution which is close to the final solution, integration 
methods like Newton's method or its variations can be used to drive the solution 
close to zero. But when the initial solution is far from the final solution, this 
method does not reach the quadratic convergence within an acceptable iteration 
count. This restriction can be partially alleviated if one uses the alternating 
directional implicit (ADI) method [76,77, 82). This method operates on the 
principle of first factoring Eq. (2.35) into easily manageable unidimensional 
operators and then sweeping once in the £ direction and next in the r\ direction 
or vice versa until a steady-state solution is achieved. If one applies the ADI 
method to Eq. (2.35) one obtains 
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[/ + AtJ^A" + <5“A + ]}aQ* = -At R h (q) 

(2.38a) 

[l + A^B'+S+J-^AQ^AQ' 

(2.38b) 


For the steady-state computation, Eqs. (2.38a) and (2.38b) are advanced by a 
local time step At which is computed for the given CFL number as [82] 


At <CFL- 


Vi + mW(I v (*h v <’'>i) 


(2.39) 


where CFL is the Courant-Friedrich-Lewy number. 

2.4 Initial and Boundary Conditions of Two-Dimensional Euler Equations 

Although the initial condition can be arbitrarily assigned, the usual practice in 
the steady-state computation is to initiate the computational domain with the 
reference values of the state variables. On the other hand, the boundary 
condition cannot be arbitrarily assigned because the solution depends on the 
unique determination of the boundary condition that stems from the physical 
reasoning of the flow. The far-field boundary condition of the Euler equations, for 
example, is derived from the propagation of information along the characteristics 
defined by the Riemann invariants as 

R ± =U±^~ (2.40) 

7-1 


Then with the help of Eq. (2.40), the normal velocities and the speed of sound at 
the body surfaces can be computed. Consequently, the inflow and respectively 
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the outflow Cartesian components of the velocities on the outer boundary can be 
computed as 

v = v. + (v-v - .«)n (2.41a) 


and 


v = v «,, + 



(2.41b) 


where 

V = [U,V] T v=[w,vf n = [Z x ,Q T (2.41c) 

In Eq. (2.41c), \ can be either in the direction or in the 77 direction. 

For inviscid flow with impermeable and adiabatic wall conditions, the 
contravariant velocity normal to the wall, and the normal derivatives of the density 
and pressure at the wall are assumed to be identically zero. Mathematically, this 
can be expressed by 

V =0 — = 0 ^ = 0 (2.42) 

dn Bn 

The extrapolated density and pressure from Eq. (2.39) are then imposed 
explicitly on the boundary and updated at every iteration until the solution 


converges. 
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2.5 Governing Equations of Quasi One-Dimensional Euler Equations 

The time-dependent, compressible, quasi one-dimensional Euler equations 
[76] are formulated both in the nonconservative and conservative form to 
elucidate two design sensitivity approaches, i.e., sensitivity equations derived 
from one-dimensional nonconservative field variables and sensitivity equations 
from one-dimensional conservative variables. 


2.5.1 Quasi One-Dimensional Euler Equations in NonConservative Form 

The quasi one-dimensional Euler equations in nonconservative forms are 
commonly expressed in terms of the conservation of mass, momentum, and 
energy. These conserved quantities of mass, momentum, and energy can be 
given, respectively, as 

^ + ^>=0 (2.43) 

dt dx 

d(puS) + 3( pS+P)S p dS =0 (244) 

dt dx dx 


and 


dt dx 


(2.45) 


For the first approach (nonconservative field variables), Eqs. (2.43) - (2.45) are 
put in their vector form as 

d(SQ) 


dt 


(2.46) 
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where 


Q = [p,pu,pe,] T (2.47) 

S is the cross-sectional area, u is the velocity, p is the pressure, e, is the total 
internal energy, p is the density, and r denotes the vector of the spatial 
derivatives. 


2.5.2 Quasi One-Dimensional Euler Equations in Conservative Form 

Unlike the nonconservative approach, this approach deals directly with the 
variation of the conservative fluxes with respect to the conservative field 
variables. To realize this objective, let us recast Eqs. (2.43) - (2.45) in a form 
amenable to conservative formulation. If one follows this procedure, the time- 
dependent quasi one-dimensional Euler equations in conservative forms are 
given as 


S(SQ) | BE 

dt dx 


-H =0 


(2.48) 


where 


H s =H s (S,Q) = 


n dS n 

0,p — ,0 

dx _ 


(2.49) 


are the source terms. 

In Sec. 2.4.1, the one-dimensional fluxes and Jacobians are computed by 
both the Beam-Warming and Van Leer flux-vector-splitting techniques. In Sec. 
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2.4.2, however, the Van Leer flux-vector-splitting techniques are only adopted. 
The discretization, and time integration of both approaches are performed as 
given in [76]. 

2.6 Initial and Boundary Conditions of Quasi One-dimensional Euler 

Equations 

For the purely supersonic flow case, the flow properties p, u. and e, are 
specified at the inflow boundary; whereas, at the outflow boundary they are 
numerically computed from the interior points using a second-order extrapolation. 
For the supersonic-inflow and subsonic-outflow conditions, the inflow boundary 
is specified as in the supersonic conditions. However, at the outflow boundary 
location, the pressure is specified whereas the density and velocity are computed 
from the interior by a second-order extrapolation. 

In the purely subsonic flow condition, the density and the velocity are 
specified at the inflow, and the pressure is specified at the outflow boundary. 
During the numerical experimentation on the various ways to specify the 
numerical and physical boundary conditions that produce stable and converged 
solutions, it was generally observed that using first-order extrapolations of the 
numerical boundary conditions both at the inflow and outflow boundaries 
produced either spurious or non-convergent solutions. However, switching the 
inflow boundary to a first-order extrapolation seemed to give comparable results. 



Chapter 3 

AERODYNAMIC DESIGN OPTIMIZATION AND SENSITIVITY ANALYSIS 


The work done with regard to general aerodynamic optimization problems is 
discussed in this chapter. The pertinent elements of aerodynamic shape design 
optimization and sensitivity analysis by the application of the variational methods 
(VM) are also presented. 


3.1 General Scope 

Since the 1950's, design and shape optimization for bodies and the numerical 
algorithms to improve analysis and optimization processes steadily grew. Now, 
with the latest supercomputers and software technology, many design problems 
are being solved in a matter of hours. Therefore, the question today is not 
whether researchers want to design numerically complex parts, but rather, which 
is the most reliable design and shape optimization tool. This quest ultimately 
culminates in devising an algorithm to get better gradient or sensitivity 
information of the functional or system responses to the change of design 
variables or control functions. 

3.1.1 Sensitivity Information 

At the present time, three main avenues exist to get the gradient information. 
The first, and the oldest one, is the finite difference approach. Depending on the 
type of approximation and number of design variables, this method requires a 
substantial number of converged solutions and is also dependent on the step 
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size, which is always arbitrary. The gradient information is commonly obtained by 
either forward, backward, or central difference approximations. With the use of 
the central difference alternative, for example, the gradient is computed as 

3d d[q(x d + /x d \x d + m d ]-d[q(x d -ax^ 

3X d ~ 2AX d 


The second category is the discrete approach. This is widely used and is 
computationally intensive. In its standard form, it also requires huge computer 
memory allocation, especially, for the two-dimensional complete Navier-Stokes 
equations and, generally, for three-dimensional equations. This can be attributed 
to the large and sparse matrices obtained from the sensitivity of the solution 
vector and grid terms with respect to the design variables. Generally, the 
sensitivity coefficient is derived [83] by the chain rule for implicit functions, for 
instance, as 


3D 

3X D 


3D } + (3D ' 


3X 


D Jq 


\ >>Q i 


ax„ 


(3.2) 


The derivative of the solution vector with respect to the vector of the design 

r)G 

variables — =- can be obtained by the implicit differentiation of descrete residual 

dX n 


which is defined as 


R(e(x J ),x J ) = o 


(3.3) 


The determinations of 


r dR ' 




and 


( 3D ^ 


\ dX * Jq 


, which are needed in the computation 


of the total gradient or sensitivitiy coefficients of the objective functional, are 
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computed by use of the chain rule for implicit functions of the discrete residual 
and objective functional with respect to the metric terms and the derivative of the 
metric terms with respect to the design variables. This process is also done for 
the gradients of the constraint functionals. 

The third way of obtaining the sensitivity gradient is through the variational 
methods. Unlike the discrete approach, variational methods derive the derivatives 
of the functional based upon the calculus of variations. The optimal control theory 
is an example of using the variational methods to derive the optimality conditions 
whereby the extremizing functions (design variables) are solved for. In general, 
variational design optimization blends the concepts of optimal control theory and 
calculus of variations [84], The optimal control theory states the conditions under 
which the control variables, parameters, and functions or the combinations of 
them can be continuously altered to meet the desired criteria. The optimal control 
theory particularly focuses on the following: (1) determination of a mathematical 
model of the dynamical or physical system, (2) determination of the admissible 
control variables or functions, (3) specification of the performance index or 
functional that can be extremized, (4) identification of the physical constraints that 
produce unique and converged solutions, and (5) construction of the augmented 
functional that consists of the objective functional and the constraints. After the 
optimal control problem is explicitly formulated, the fundamental principles of 
calculus of variations determines the variation of the augmented functional to the 
variation of the admissible control variables or functions or design variables. 

3.1.2 Discrete Versus Variational Methods 

The discrete approach describes the local behavior of the functional that is to 
be differentiated in an infinitesimal region. The variational approach, on the 
otherhand, integrates the functional that is to be differentiated over the 
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continuous domain or boundary. While both procedures are mathematically 
equivalent, the discrete approach is based on differentiation and the variational 
approach on the integration. Conceptually, differentiation is an inductive process, 
whereas integration is a deductive process. Computationally, integration is a 
smoothing operation, i.e., many weak integral (variational) forms can be stated, 
whereas differential (discrete) approach is a noisy operation. In short, the 
discrete sensitivity analysis involves prior discretization of the continuous flow 
field equations and boundary conditions before they are differentiated with 
respect to the design variables. On the other hand, in the variational sensitivity 
analysis, the continuous state equation in the weak form of the integral is first 
differentiated by use of the principle of calculus of variations and optimal control 
theory to derive the costate equations, transversality conditions, and functional 
sensitivity derivatives, and then they are discretized. 

3.2 Constrained Optimization Methodology 

A constrained optimization method in general encompasses three elements of 
optimization, i.e., design variables, constraints, and objective function. 

3.2.1 Design Variables in Variational Sense 

In most aerodynamic optimization problems, the design variables are 
generally of a geometric nature, such as the coefficients of some geometric 
functions, surface grid points [83], aerofunctions [85], or polynomial functions 
such as Bezier-Bernstein functions [42, 86] and spline functions [87], 

Variational methods treat the boundary of the domain in a continuous fashion, 
and therefore, the boundary is considered as part of the solution to the design 
problem. With the assumption that the domain ft is sufficiently regular, the 
location of points on the boundary X r can be considered as a continuous design 
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variable (Figs. 3.1 and 3.2). Mathematically, the coordinates of the varying 
boundary in the continuous sense can be expressed as 

X r =f(X D ) (3.4) 

where X d are the design variables. In aerodynamic optimization problems, the 
vector of design variables is provided for very limited and simplified geometries, 
for instance, 4 digit NACA airfoils and some nozzles. However, for general- 
purpose geometries, these control points must be determined through iterative 
methods from certain functional relationships such as the Bezier-Bernstein 
polynomials [42], Because these polynomial functions are known to generate 
smooth curves and surfaces for a minimal number of control points, the function 
/ which describes the curve for the two-dimensional problem, is given by [86] 

f(u)=ix di B IH (u) for u e [0,1] (3.5) 

1=0 

where 


B in (u) = C(n,i)u‘(\-ur > 


(3.6) 


C(n,i) = 


n\ 

/!(«-/)! 


(3.7) 


In Eqs. (3.5) - (3.7), B in (u ) are the blending functions, which are key to the 
behavior of the curve, C(n,i) are the binomial coefficients, u is the normalized 
arc length and n is the order of the Bezier-Bernsten polynomials. In this study, 
with the use of Eqs. (3.4) - (3.7), the location of the control points can be 
considred as the design variables. 
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3.2.2 Constraints 

Constraints are the integral parts of the optimization procedure that influence 
the final outcome of the functional. They can be geometrical, flow-type, equality 
or inequality constraints, or a combination of all or some that depends on the 
particular optimization problem one wants to address. 

In the design optimization process, certain constraints are bound to be 
satisfied while the others are violated. Those which are satisfied encompass the 
feasible domain, while the violated constraints belong to the infeasible domain. 

In the variational formulation of design optimization problems, the flow-type 
constraints are expressed in the integral forms. The geometrical and side 
constraints, on the other hand, can be formulated either in the integral or discrete 
forms. For the general variational approach, generic flow-type constraints are 
expressed as 

G j (Q,n) = jg j (Q,n)dT<0 for; = 1 . 2, ..., nconf (3.8) 

f 

where F, is the deformed boundary and nconf is the number of generic fluid-type 
constraints. The generic geometric-type and the side constraints can also be 
given as 

Gj(X d ) < 0 for j = nconf+1 , nconf+2, ..., neon (3.9) 

and 

X'™‘ r < X iD < X^ r for / =1,2 ndv (3. 1 0) 
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where neon is the total number of constraints, and ndv is the number of design 
variables, respectively. 

3.2.3 Objective Functional 

In the variational methods (VM), the objective functional is defined in the form 
of a definite integral involving an unknown state function Q, which can be 
dependent on some normal vectors n and other problem parameters. Then, the 
objective functional is extremized at the converged state solution over the curve 
of the surface described by the vector of design variables. Mathematically, a 
generic functional on the boundary J f , is defined as 

J T (Q,n) = \D(Q,n)dT ( 3 . 11 ) 

f 

where D, for the two-dimensional problem, is the objective function specified on 
the curve or boundary. The selection of the objective function is mostly dictated 
by the flow physics. 

3.3 Variational Formulation of Aerodynamic Optimization Problem 

In a design optimization where the constraints are absent, the necessary 
condition for the objective function to reach its optimal solution is when the partial 
derivatives of the function with respect to the design variables are all identically 
zero. The sufficient condition for optimality can be further augmented by requiring 
the Hessian matrix of the function to be positive-definite at the design points. 
Realistic optimization problems involve constraints which have functional 
relationships to the functional through the design variables or state variables. 
When constraints are involved in the optimization problem, the partial derivatives 
of the functional and the constraints cannot be zero at the same time since they 
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are functionally related to each other through the optimality criteria [84, 88]. One 
common practice is to cast the constrained optimization to unconstrained 
optimization through the introduction of the weighting functions or Lagrange 
multipliers X{X). The other is to sequentially solve a linear or quadratic 
programming problem, which is an approximation of the original constained 
minimization problem. In the later approach, one needs to derive the sensitivities 
of the performance functional. 

In the following paragraph, the objective functional defined in Eq. (3.11) will 
be used as an example to facilitate the discussion of the procedure to derive the 
aerodynamic sensitivity equations by the variational methods. To start the 
derivation, the steady state solution of Eq. (2.6), i.e., the residual R(Q), is written 
as 


R(Q) = E X +F y = 6 (3.12) 

and the generic boundary conditions are expressed as 

H{Q,n) = 0 (3.13) 

Without changing its value, the objective functional J f can now be modified as 
J T ^Q,nXii) = J f +\^{XmQ)dO. + \p T H(Q,R)dT (3.14) 

ti r 


where r and Q are the deformed boundary and domain, and X and p are 
vectors to be determined. 
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3.3.1 Standard Formulation of an Aerodynamic Optimization Problem 

A mathematical formulation of the constrained optimization problem can be 
expressed as 

min {/ r } (3.15) 

subject to 

Gj(an) = j gj (Q,n)dr< 0 for/ = 1 , 2, nconf (3.16) 

r 

Gj(X d )< 0 j = riconf+1, nconf+2, ..., neon (3.17) 

and 

X£" < X tD < X t T” / =1.2, .... NDV (3. 1 8) 

where the flow field variables Q are the solution to the state equations, R(Q). 

The problem statement clearly indicates that the state equations are the 
integral parts of the optimization process and, therefore, must be represented by 
the highest level of flow field approximations and solved by the most efficient 
numerical techniques as dicussed in Chap. 2. 

3.3.2 Derivation of Functional Sensitivity Equations 

In the derivation of the sensitivity derivatives for the functional and 
constraints, the spirit of Ref. 65 is still kept in implementing the variational 
methods . But unlike Ref. 65 the variation is performed on the conservative 
variables and fluxes with respect to the variation of a variable domain. Let us 
define the following expressions for later use: 
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dE_ 

dQ 


= A 


and 



(3.19) 


to be the Jacobian matrices in the x and y directions. The variation of the fluxes 
(to the first order) can be written as 

SE =ASQ and 8F=B5Q (3.20) 

Then the fluxes on the deformed space due to the variation of the boundary can 
be approximated as 

E(U) = E(Q) + 8E(Q) and F(U ) = F(Q) + SF(Q) (3.21) 


By application of the principles of calculus of variations and use of the results in 
Appendices A and B, the variation of the modified functional can be 
approximated by [84] 

= ( 3 - 22 ) 


where X and X are position vectors of the deformed and undeformed coordinate 
systems, respectively. Then, with Eqs. (3.19) - (3.21) and Eqs. (A.2), (A. 6), and 
(B.5) from Appendices A and B, the Taylor expansion of the integrand of Eq. 
(3.22) [84] is computed (the linear part relative to e ) as [84,88,89] 
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8J ra = J \d + D q 5Q + D n 8n + D Kdn^dT 

+ j(fi T + 5/2 T ){/7 + H q 8U + H n 8n + //Ktf/jJdT 

r 

+ J |(X T + 8X")\E x {Q) + (8E(Q)) X + F y (Q) + (<5F(0),]}|-/s|^ 
- |A r [^(0 + ^(0)]^-J[^ + /i r ^]jr 


(3.23) 


where e is a small parameter, J s is the space transformations matrix that is given 


in Eq. (B.5) as |7 5 | = 


I + V.8X 


and the quantity k in Eq. (3.27) is the curvature 


and can be calculated as [88] 


k = -V o n 


(3.24) 


where ° denotess a dot or inner product and n is the unit normal which can be 
computed from the grid generating routine or from the analytical derivatives of the 
Bezier-Bernstein polynomials as [86] 




3 " 

du du 1 

du 



(3.25) 


where <8> is a vector multiplication sign and f{u) is defined in Eq. (3.9). In Eq. 
(3.23), 8n is defined as 

8n = 8X°n ■ (3.26) 


Now, by taking only the linear terms of Eq. (3.23), one obtains 
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8J ra = j \d q 5Q +D H 8n + DK8n + fi T ^H Q 5Q+H n 5n + HKdn^dT 

+ \{l T [(5£(Q)), + (5F(G)),]}dQ + J |(5A) X [£ x «2) + F,(Q)]}dQ 

o n ^ J 

+ j8fi T Hdr + J X T [E x + F y ]8ndT (3.27) 

r r 


With Eq.(3.20) and performing integration by parts, the second term in Eq. (3.27) 
is expressed as 


\{X T [(8E(Q)) x + (8F(Q)) y ]}dCl= \{-X T x A-X T y B)8QdSl + 
n a 

\{£ T [An x +Bn y ]S§\dr (3.28) 


Where A and B are defined in Eq. (3.19). Substitution of Eq. (3.28) into Eq. 
(3.27) gives 

8J Ta = \\D Q 8U + D n 8n + DK8n + fi T \H Q 8t) + H n 8n + HKSt^d^ 

+ \{i T [An x +Bn y ]8t)\ir + \[-X T x A-X T y B]8UdQ. 

r n 

+ j{( 5 *ffc(2 ) + ^(G)]W+ \Sfi T Hdr+ \jj[E x + F y )8ndr (3.29) 

n *■ J r r 


Note that for the arbitrary variations of 8X and 8jd and with Eqs. (3.12) and 
(3.13), the last two terms in Eq. (3.29) are identically zero. Then Eq. (3.29) 
reduces to 
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8J Ta = \^D q 8Q + D n 5n + DK8n + fi T ^H Q 8Q + H n 8n + //*r5«]|rfr 

+ \{x T \An x +Bn y ]^^+\[-X\A-X\B^8Ud^ + \X T [E x + F y ]8ndT 
r 1 ci r 

(3.30) 

In Eq. (3.30), the vectors A, and p can now be determined to eliminate the 
terms associated with 8Q . Consequently, the costate (adjoint) equations are 
given as 


-A T X x -X y B T = 0 in Q (3.31) 

Upon the combination of Eqs. (3.30) and (3.31), the variation of the functional 
becomes 

8J ra = J j D q 8U + D n 8n + DK8n + P t [h q 8@ + H n 8n + //x-foJJdT 

+ jlX T [An x +Bn y ]8U]dT + jX T [E x +F y ]8ndT (3.32) 

r 

With Eq. (A.1 1) in Appendix A, we now express 8Q in terms of SQ to get 
8U = 8Q-?S-8X (3.33) 


For the sake of computational simplification, the variation of 8X on the boundary 
is limited only to the y component in this study, i.e., 


a=[0,s>f 


(3.34) 
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and Eq. (3.33) is simplified as 

S§ = 8Q~^8y (3.35) 

dy 

Also an approximatation of Eq. (3.26) and use of Eq. (3.34), Eq. (3.26) can be 
written as 


8n = 8X o n 


= [0,<5yf °[n x ,/i y ] 


= n y °8y 


(3.36) 


By use of Eqs. (3.33) - (3.36), Eq. (3.32) is now given as 

8J ra = J {D Q SQ + n y D„8y + n y D >c8y - n y D Q Q y 8y +}dT 
r 

+ J {[i T [H Q 8Q + n y H n 8y + n y HK8y - n y H Q Q y 8y]]dr 

r 

+ ^}J^An x + Bn y Y8Q- ti y Q y 8y)^dr + J X t ^E x + F y ^i y 8ydT (3.37) 
r r 

For arbitrary 8Q and the variation of y on the boundary r, Eq. (3.37) gives the 
boundary conditions for the costate equations and the sensitivity equations, 
respectively, as 

^An x + Bn y f X+Dl+H^ = 0 on T (3.38) 


and 
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S/r. = J{i>. + Dk- D Q Q y +fi r [H, + HK-H/^noydT 

r 

- \{X T [An x +Bn y ]Q y }n y 8ydr + \^ t \E x +F^i y SydT (3.39) 

r r 


The unique determination of Eq. (3.39), therefore, demands the unique and 
converged solutions of Eqs.(3.12) and (3.13), (3.31), and (3.38). 

3.3.3 Derivation of Constraint Sensitivity Equations 

With the constraints defined in Eq. (3.16), the residual, and the boundary 
conditions, Eqs. (3.12) and (3.13), one can formulate the modified constraints as 

°>r.(a *. ■ hSt ) ■ = I */(< 2.«K + + J ji/ff (asK 
r n r 

for j =1,2 nconf (3.40) 

By following the same procedure as was done for the objective functional, the 
costate equations, boundary conditions, and the constraint derivative coefficients, 
respectively, can be expressed as 

-2PX^j = 0 for j= 1, 2, ..., nconf in Q (3.41) 

{[Jrt, + +4 + = 0 for y = 1.2 nconf on r (3.42) 

SGj r , = J {*>. +«A"*aA + ''/[«. + 

r 

- + + JiJ[F,+F,]ti,$yrfT 

r 


for j =1,2 nconf 


(3.43) 
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As can be discerned from Eq. (3.43), the computation of the constraint sensitivity 
equations requires the solution of a new set of costate equations and boundary 
conditions as many times as the number of constraints. 

3.4 Numerical Optimization 

In Sec. 3.2, the elements of gradient-based constrained optimization, i.e., the 
parameterization of the boundary, objective and constraint functionals, and the 
sensitivity derivatives for an aerodynamic design optimization problem, are 
presented. The next logical step is to specify a numerical optimization technique 
to search for a better design. The feasible direction method developed by 
Vanderplaats and Moses [90] and used by Haftka and Gurdal [91] is 
implemented in our study. Two steps are essentially followed in this approach. 
The first step is to determine the search direction, S , and the second is to 
compute the magnitude of the step size a. These two quantities can be 
computed as proposed in Refs. 90 and 91 A typical computation of the feasible 
direction starts at the boundary of the feasible domain, and its magnitude and 
directions are kept constant as long as the search direction keeps the design 
variables in the feasible domain while improving the performance index. 
Otherwise, a new search direction and step size are recomputed with the new 
gradient information and this process continues until the optimality is met. 
Mathematically, the feasible direction can be formulated as 

S T oVg.<0, (3.44) 

where / is part of the active constraints and the usable direction at a point is 
given by 



The change in design must be sought along the combination of the useable and 
feasible directions so that the functional or the performance index is reduced as 
much as possible, and the design is kept away from the constrained boundary as 
much as possible. By use Eqs. (3.44) and (3.45) in the method of feasible 
directions, the new design variables are updated as 

Xp + 1 = + aS (3.46) 

where n is the iteration number. The values of the design variables are 
continuously altered until the criteria for the optimal solution of the performance 
index are satisfied. 
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Fig. 3.1 Variation of domain by a One-parameter family of mapping 


4 


fi 


20 



Fig. 3.2 Boundary variation normal to the original boundary, 20 
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Chapter 4 

COSTATE EQUATIONS AND SOLUTION METHODS 

4.1 Introduction to the Numerical Integration of Costate Equations 

The coefficients of the costate equations are constant matrices whose 
components are derived from the converged solution of the state equations. (See 
Appendix E). They are globally constant in time and locally constant in space. But 
the interpretation of constant matrices must be understood in a sense that, 
during the time integration of the costate equations, only the costate variables 
evolve in space and time to convergence. The costate equations are identical to 
the Euler equations in form, but mathematically, they are different in the sense 
that they do not meet the homogeneity requirement to put them in a conservative 
form like the Euler equations. From the numerical view point, one can adopt any 
solution algorithm, which is used for the Euler equations, to the costate 
equations. This can be explained by the fact that the fluxes on the cell faces or at 
grid points can be artificially constructed by approximating the solution vector of 
the costate equations either on the cell face from the right and left sides of the 
cell centers or at the grid points in exactly the same way one does for the state 
fluxes and solution vector. 

4.2 Costate Equations 

The costate equations, like the state equations, are solved by use of the time 
dependent techniques. The Eqs. (3.31) and (3.41) are, therefore, modified to 
include the unsteady term with the proper signs so that this time dependent 
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technique is fully exploited. Thus, for instance, Eq. (3.35) in the generalized 
coordinates system is expressed as 

±\ t -A t X 4 -B T X n =6 (4.1) 

The proper sign selection of the time term is dependent on the complementary 
property of the well-posed boundary conditions of the state and costate 
equations. For Eq. (4.1) to be well-posed, the positive sign of the time term is 
selected, and Eq. (4.1) becomes 

(4.2) 

4.3 Boundary (Transversality) Conditions 

The objective functional boundary conditions, i.e. Eq. (3.38), in their general 
forms are again for the sake of convenience presented here as 

{< 8Q T [An x + Bn y ] T X + D t q +H^ = 0 on T (4.3) 

The objective functional and the no-mass penetration conditions are defined only 
on the solid boundary, and hence their derivative contributions in Eq. (4.3) are 
identically zero. Therefore, the boundary conditions for the inlet, exit, and center- 
line reduce to 

^8Q T [An x +Bn y ] T ^ = 0 on (r^r^r.*) (4.4) 

For the supersonic flow flow, the inlet condition is known, and hence the 
variation of the vector of the flow field is identically zero. Therefore, with Eq. (4.4), 
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the values of the costate variables at the inlet boundary can be approximated 
from the internal stencils. Because the the vector of the flow field is computed 
from the internal grids at the exit plane, Eq. (4.4) gives four linear independent 
equations for the costate variables, which result in all the costate varibles to be 
identically zero. On the centerline, the normal velocity is known to be zero, and 
one of the costate variables, for instance A 3 , is assigned a value, and the 

remaining flow field quantities are to be detemined from the resulting 3x3 
system of equations as given in Eq. (4.4) 

One way of treating the boundary conditions, i.e. Eq. (4.3), is to use Eq. (3.13) 
and to find a relationship between the conservative field variables Q by taking the 
variation of Eq. (3.13). This procedure eliminates the constant Lagrange 
multipliers p and modifies the functional sensitivity derivatives, Eqs. (3.39) and 
(3.43) by a term resulting from the variation of the normal vector h at the solid 
boundary. 

On the solid boundary, on the other hand, the costate variables are 
determined by use of the complete form of the compatibility relationships and the 
sign of the eigenvalues of the costate Jacobian matrices. Once the values of the 
costate variables on the solid boundary are computed, the constant Lagrange 
multipliers p of the no-mass penetration condition can be calculated by solving 
the complete set of the boundary condition. The results presented in this study 
are obtained by solution of the complete boundary conditions as given in Eq. 
(3.42). 


4.4 Linearization of Costate (Adjoint) Equations 

By the same linearization procedure we used for the state equations, Eq. (4.2) 


can be approximated as follows: 
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X, -At-^-{a t X ( + B 7 i„j = {.4 r A, + B’X„y 

(4.5) 

X, - Ar£{[.4 r S f + B^„]X} = {X 7 X { + S 7 X„}' 

(4.6) 

v - + - {**« + 

(4.7) 


By approximation of the time and space terms, Eq. (4.7) then becomes 
AA" - At [a t <^ + 5 T <5„]AA" = At{[a t ^ + (4.8) 

4.5 Time-Integration Methods 

In this study we have used both the implicit and explicit, i.e., the ADI and 
Runge-Kutta time-integration, respectively, methods to drive costate equations to 
steady state. For the implicit method, the ADI factorization of Eq. (4.8) is used to 
split it into the £ and r\ sweeps. Let us define the right side R x of Eq. (4.8) as 

R x = At{[(A-) t ^ + +{A + ) T 8- +(ET) t 5; + (S + ) 7 '5;]a}" (4.9) 

where R x is the residual for the costate equations. Also, Eq. (4.8) can be put in 
its split form of Jacobians and fluxes as 

{/ - Ar[( A+ ) T 8~ + ( AT ) T 8\ + (B + ) T 8~ v + (B~ f 8 + v ]}aF = R x (4.10) 


Then the E, and rj sweeps of Eq. (4.10) are given as 
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{/ - Ar[(/4 + ) r 5' + (/T) r <5j]}AA n ‘ = R x (4.1 1 ) 


and 

{/ - Ar[(B + ) t <5; + (B~ ) T 6 ; ]|aA" = AA n * (4.12) 

For the explicit method, the four stage Runge-Kutta algorithm [81] is adopted 
to compute the vector of the costate variables as 


II 

(4.13) 

A (2) = A" + a 2 AtR x (1) 

(4.14) 

A (3) = A" + a 3 A tR x (2) 

(4.15) 

A <4) - A" + cc 4 AtR x {3) 

(4.16) 

A" +1 = A" + + 2/? i 2) + 1R x ] + R x } ] 

(4.17) 

where 


a 2 =y 2 « 3 = l /i <* 4 =1 

(4.18) 


In advancing the costate equations in time, both the implicit ADI and explicit 
Runge-Kutta methods are employed. The results, however, will be presented 
only for the ADI method. 
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Chapter 5 

STABILITY ANALYSIS OF THE TWO DIMENSIONAL COSTATE AND EULER 

EQUATIONS 

5.1 Rationale for Stability Analysis 

The most popular schemes to advance the Euler equations to steady-state 
solutions are, among others, the implicitly factored time-integration schemes. 
Approximate factorizations unfortunately introduce errors which propagate 
throughout the computational domain. As a result of this, the stability limit is 
drastically reduced, and the convergence rate deteriorates. To propose the range 
of the CFL numbers for which the allowable maximum eigenvalues are predicted, 
a stability analysis of the Euler and costate equations are conducted. 

5.2 Introduction to Stability Analysis 

In solving coupled time-dependent partial differential equations like the Euler 
equations, one takes advantage of the hyperbolic nature of the equations. 
Because of this fact, numerical upwind methods are devised according to the 
direction of the flow information along the characteristics. The most common 
upwind methods that take into account the hyperbolicity of the equations are the 
Beam-Warming, Van Leer vector-splitting, and the Roe flux-differencing methods 
[92], On the basis of the eigenvalue decompositions, the fluxes and Jacobians of 
these methods are split into the backward and forward contributions. The Van 
Leer vector-splitting method [80] is shown to produce sharper shocks than the 
Beam -Warming flux - vector - splitting approach [79]. In this analysis, the 
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Beam-Warming splitting for the one-dimensional Euler equations has been used 
for the purpose of comparison only. Otherwise, the Van Leer vector-splitting 
technique has been adopted throughout for both the one-dimensional and two- 
dimensional Euler equations. 

As mentioned in Sec. 1.2, Jesperson and Pulliam [71] studied the stability 
characteristics of the Euler equations for different flux-splitting methods. 
Anderson and Thomas [72] further conducted stability analysis on the complete 
three-dimensional Euler equations. In their analysis, they have investigated 
specifically three kinds of splittings: three-factor spatially split, two-factor 
eigenvalue split, and two-factor combination split. All three splittings have 
different levels of factorization errors. During derivationof the generalized 
complex eigenvalue equations, they have also used first-order differencing on the 
implicit side (leftside) and fully upwind second-order differencing for the residual 
(right side) part. For a Mach number of 0.8 and an angle of attack of 0°, they 
found out that the three-factor splitting has lower CFL stability (CFL = 20), 
whereas the other two-factor splittings are stable for all CFL numbers considered. 

In their pursuit to optimize the PROTEUS code with multigrid methods, 
Demuren and Ibraheem [73] have conducted an extensive and complete stability 
analysis of one-, two-, and three-dimensional Eluer and two-, and three- 
dimensional Navier-Stokes equations . They have looked at not only the ADI 
factorization but also the LU approximation factorization for Euler equations and 
Navier-Stokes equations with various levels of dissipation terms. These 
inclusions, in fact, encompass the most recent and commonly used 
approximation numerical methods, specifically, the upwind and central difference 
approximations. Their conclusions are in line with Ref. 72, i.e., the CFL range 
over which the maximum eigenvalues are optimized decrease as the 
dimensionality of the problem increases. The stability deterioration is dependent 
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or augmented by the type of discretizations and factorizations employed in the 
numerical computations. 


5.3 Objective and Motivation 

In the preceding section, the main features of different schemes are analyzed 
with the help of the amplification factors for the complex eigenvalue boundary 
equations [71 - 73], The objective of the present stability analysis is not to revisit 
the complete stability analysis of different upwind schemes of the Euler and 
costate equations with varying approximations. Rather, the main thrust is to 
investigate the stability analysis of the costate equations with only the spatial 
upwind factorization scheme and the Euler equations are included for 
comparison reasons. This procedure is necessary because the costate equations 
are similar to the Euler equations, and consequently, one will adopt the same 
numerical technique for the costate equations. The numerical stability and quick 
convergence of the costate equations are very detrimental because the 
computation of the sensitivity gradient is directly influenced by the converged 
solution of the costate equations, which ultimately controls the whole analysis 
and optimization process. For this obvious reason, one has to investigate the 
CFL range over which the maximum, average, and L2-normed eigenvalues are 
extremized. 

Like the previous researchers [71 - 73], first-order difference approximation on 
the left side and second-order upwind differencing on the right side of the Euler 
and costate equations are used. To investigate the stability limits of the PDE for 
the stated factorization of the Van Leer schemes, one solves the generalized 
complex eigenvalue boundary value problems given in Secs. 5.4.1 and 5.4.2. In 
doing so, one computes the maximum, average, and L2-normed eigenvalues in 
the range of 0<co x ,co y <2 n for the desired series of CFL number. Also the 
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smoothing factor cr, which is related to the damping of the high frequencies in the 
multigrid methods, is computed as the absolute value of the maximum 
amplification in the range of -^ < max(a Following Anderson and 

Thomas [72], the time step which is used in the computation of the amplification 
factors for the two-dimensional case is given by 


A t<CFL\ 



- 1-1 


1 


1 


Ax 2 Ay 2 


(5.1) 


5.4 One- and Two-Dimensional Euler Equations 

The one- and two-dimensional Euler equations in conservative forms and in 
Cartesian coordinate systems are given first. Then, they are discretized in the 
upwind fashion depending on the positive and negative fluxes and Jacobians. For 
the one-dimensional Euler flow, the source terms are included to investigate their 
eventual influence in the stability limit of the flow characteristics. Finally, the 
approximations and discretizations of the one- and two-dimensional costate 
equations in Cartesian coordinate systems are given for completeness. 


5.4.1 One-Dimensional Euler Equations 


f(fg) | BE 

dt dx 


-H= 0 


(5.2) 


By an application of the Euler backward approximation in time on Eqs. (5.2), one 
obtains 


5/ + Ar • 


dx 


(A' + A-)-B h 


j&2 = -A; ■ 




(5.3) 
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{5/ + Ar ■ [(8;A + + 8 + x A')- B h ]]AQ = -At • [( 8~ X E + + 8 + x E ~ ) - H s ] (5.4) 

The use of E ± = A ± Q \ into Eq. (5.4) results in 

{5/ + Ar • [(8; A* + 8 + x A-)~ B h ]]aQ = -Ar • [<S;(a + <2) + 8 + x (a~q) - H J (5.5) 

where the fluxes, the Jacobians, and the source terms are given in Appendix C. 


5.4.2 One-Dimensional Costate Equations 

From Eq. (2.48), the one-dimensional costate equations in Cartesian 
coordinate systems are derived ( see the detail in Chap. 6) as 

X, - A T X x - BjX + SC T = 6 (5.6) 


where B h , and C are given in Appendix C. By the same procedure as for the 
Euler equations (Eqs. (5.3) - 5.5)), the following approximations for the costate 
equations are obtained. 


{/- Ar[5;(A + ) r + 5 I + (A-) T ]}Ar = Ar[<5;(A + ) T (5.7) 


Although the Jacobian matrices of the costate equations are transposed, they are 
given in Appendix C. 

5.4.3 Two-Dimensional Euler Equations 

From Eq. (2.7), the conservative two-dimensional Euler equations in 
Cartesian coordinates are expressed as 
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dt dx dy 


(5.8) 


By application of the Euler backward approximation in time on Eqs. (5.2) and 
(5.7), one gets 


/ + Ar • 


-?-(A + +A-) + ^-(B + +B-) 
dx ' dy' ' 


A Q n = -At ■ R n 


(5-9) 


where 

R n =4-( E+ + E ~) + -ir( F+ + r ) (5.10) 

dx dy 

Then, further representation of Eqs. (5.4) and (5.5) in the delta form gives 
{/ + Ar ■ [(5;A + + 8* x A ' ) + (8;B + + <5;fl-)]}A(2" = -Ar • R n (5. 1 1 ) 

where 

R H = (8~E + + 8+E-) + (5;F + + 8;f~) (5.1 2) 

By use of E ± = A ± Q\ and F ± =B ± Q in Eq. (5.12), Eq. (5.12) assumes the form of 

{/ + Ar • [8;A + + 8*A~ + <5;B + + <5;JT ]}aQ" = 

-Ar - {<5;(A + <2) + 8:{A-Q) + 8;{B*Q)+8;{B-Q)) n 


(5.13) 
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Factorization of the left side of Eq. (5.13) results in: 

{/ + A! + <SX)}{; + At (O.fi* + s;b - )} AQ- = 

-At ■ {s;( A*Q)+ s:(a-q) + 5;(b’q) + <5;(S-e)}" 

(5.14) 

5.4.4 Two-Dimensional Costate Equations 

With Eq. (3.31) and addition of time term with its appropriate sign, the two- 
dimensional costate equations are 

X,-A T X x -B T X y =0 (5.15) 

Adoption of the same method as for the two-dimensional Euler equations (Eqs. 
(5.9) and (5.13)) results in 

{/ - Ar • [<5;(A + f + <5;(/T) T + 8;(b + ) T + ^(zrfJjAA" = Ar • R n x (5.1 6) 

where 

K = {,s;[(M*)^] + 5;[(/i-) r i] + 5;[(B*) r A] + «;[(B-fA]}" (5.17) 

Similarly as for Eq. (5.13), factorization of the left side of Eq. (5.16) gives 

+ ^*[(4-)'a] + 5;[(B*) r a] + b;[(b-) t i]}' 


(5.18) 
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For the time increment of the flow and costate field variables, respectively, 
the following approximations are used: 


AQ n = Q n+] -Q n 

(5.19) 

Ar = A" +1 -A" 

(5.20) 

5.5 Solution Algorithm 

To perform the Von Neumann stability analysis of the Euler and costate 
equations, the flow field quantities are considered to be constant in the Cartesian 
coordinate system. For a constant flow field, the Jacobians are also assumed to 
be constant. With this assumption, the flow quantities can be expressed in terms 
of wave fronts [72, 77], Thus, the wave fronts for the Euler and costate flow 
variables are represented, respectively, by 

Q n = v n P 0 e H “°'* Ja>,) 

(5.21) 

and 


r = <p n P 0 e ni<o,+it>,) 

(5.22) 


where / is the imaginary number defined by / = V- T and co x , co y are the x and 
y modes, and <p , P 0 are the amplification factors and initial constant vectors, 
respectively. Now Eq. (5.21) is substituted into Eqs. (5.5) and (5.14), and Eq. 
(5.22) into Eqs. (5.7) and (5.18). If the resulting equations are simplified, then the 
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generalized eigenvalue equations of the Euler and costate equations are 
obtained, and these are given in the next subsections. 

5.5.1 Complex Eigenvalue Equations of Euler and Costate Equations 

The one- and two-dimensional complex eigenvalue boundary problems of the 
Euler and costate equations can be generally represented as 

[L±k)v = L(pV (5.23) 


where L and R are the left and right side Fourier symbols and V are the 
eigenvectors. Here the positive sign corresponds to the costate equations and 
the negative sign to the Euler equations. The various matrices defined in the one- 
dimensional Euler and costate equations are given in Appendix C, and the 
matrices for the two-dimensional analysis are given in Appendix D. 


5.5.2 One-Dimensional Euier-Fourier Symbols 

The insertion of the one-dimensional form of Eq. (5.21) into Eq. (5.5) gives the 
one-dimensional Fourier symbols of the Euler equations as 

L = SI — At B h + j^-[(A + “ )(1 “ cos co x ) + (A + + A~)l sin « x ]| (5.24) 


R = - ^ )(3 + cos2ft) x - 4cosft) x ) + (A + + A~)(4sin<y x - sin 

-AtB h 



(5.25) 
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5.5.3 One-Dimensional Costate Fourier Symbols 

The insertion of the one-dimensional form of Eq. (5.22) into Eq. (5.7) gives the 
one-dimensional Fourier symbols of the costate equations as 

L = I - {^“[(^ " ^") T (l “ cos + (' 4+ + '4') r / sin ..]} (5.26) 

R = -A - ) 7 (3 + cos 2 ft)* -4cosry i ) + (/t + + /T) r (4sin co x - sin2ft)J/J 

+ A tBl (5.27) 


5.5.4 Two-Dimensional Euler-Fourier Symbols 

The insertion of Eq. (5.21) into Eq. (5.14) gives the two-dimensional Fourier 
symbols of the Euler equations as 



(5.28) 


R = ^-[(i4 + -/l")(3 + cos2c(; :t -4cosco x ) + (/^ + + A )(4sino> x -sin2<y x )/] 

+ - B “)(3 + cos2c0 > -4cos£y y ) + (B + + B~)(4sinG) y -sin2a) y )/] 


(5.29) 


5.5.5 Two-Dimensional Costate Fourier Symbols 

The insertion of Eq. (5.22) into Eq. (5.18) gives the two-dimensional Fourier 
symbols of the costate equations as 
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L = {l _ ^[( A+ -A~) r (l-cos<yJ + (/t + + A') r /sin<w i ]| 
jl - [(fi + - B~ ) T (l - cos a y ) + (B + + B~ f / sin co y j 


(5.30) 


) r (3 + cos2o> x -4cosa) x ) + (A + + A ) (4sin co x — sin2o> x )/j 
^ J(£ + -B~) T {3 + cos2co y -4coso) > ) + (fi + + B~f (4sin Ct> > -sir^cyj/j 

(5.31) 


2Ay 


5.6 Results, Discussion, and Recommendations 

To confirm that factorization indeed puts a restriction on how to choose the 
CFL range for the stable solution of the approximated PDE, a stability 
investigation of two-dimensional Euler equations in unfactored form is performed. 
Figure 5.1 clearly depicts that the maximum allowable amplification factor is 
stable at all CFL numbers considered. As the case of solving Euler equations in 
unfactored form is computationally prohibitive, one would rather revert to solving 
the factored and easily invertible operators in a unidimensional mode. 

To gain insight and confidence in solving numerically adjoint equations for 
optimization purposes, a systematic stability analysis is performed on the one- 
and two-dimensional adjoint equations along with the corresponding one- and 
two-dimensional Euler equations. In addition to the usual analysis of the quasi 
one-dimensional Euler equations without the source terms, stability 
characteristics with the source terms are also conducted. The study shows that 
the range of the CFL number for which the upwind schemes are stable is not 
basically affected by the source terms (Fig. 5.2). Note that the costate equations 
have two distinguishing features. The first one is that the transposed Jacobians 
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are constructed from the converged solution of the Euler equations. This 
indicates that those matrices are globally constant in time even though they vary 
in space. The second feature is that the matrices have negative entries. 
Mathematically, and from the stability view point of a matrix, the eigenvalues of a 
transposed matrix are no different from the eigenvalues of the nontransposed 
matrices. Therefore, transposition of the matrices has no bearing on the stability 
characteristics of the costate equations. Rather, the deterioration of the stability 
limit of the costate equations stems from the negative sign of the transposed 
matrices (Figs 5.4 and 5.5). Therefore, the correct direction of the flow for the 
costate equations is imperative. As can be discerned from Figs. 5.2, 5.3, 5.6, and 
5.7, the trends of the stability of the adjoint equations are similar to those of Euler 
equations if the flow field is transformed as presented in Ibrahim and Baysal [65]. 
When the direction of the flow is not taken into account, which is the case with 
the costate equations, the stability is limited to very small CFL numbers (Figs. 5.4 
and 5.5). This low-stability problem becomes especially apparent when there is 
high flow discontinuity in the flow field as mentioned in Ibrahim and Baysal [65]. 
Therefore, a transformation of the flow field or another numerical approach to 
account for the reverse flow direction of the costate equations is highly 
recommended. 
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1 0 1 CFL (log) 2 

Fig. 5.1 Amplification of unfactored 2-D Euler equations. 


Amplification 


64 



-1 0 i CFL (log) 2 3 


Fig 5.2 Amplification factor for 1 -D Euler equations. 
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-10 12 3 

CFL (log) 

Fig. 5.3 Amplification of factored 2-D Euler equations. 
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Maximum allowable (=1.0) 
Smoothing factor 
L2-normed 
Maximum 


O %p 


- 3 - 2-10 1 2 3 

CFL (log) 

Fig. 5.4 Amplification of factored 1 -D costate equations without space transformation. 
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-10 12 3 

CFL (log) 

Fig. 5.5 Amplification of factored 2-D costate equations without space transformation. 
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-10 12 3 

CFL (log) 

Fig. 5.7 Amplification of factored 2-D costate equations with space transformations. 
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Chapter 6 

DESIGN OPTIMIZATION OF INTERNAL FLOWS USING QUASI 
ONE-DIMENSIONAL EULER EQUATIONS 

6.1 Introduction to One-Dimensional Design Optimization 

Two approaches, one based on the nonconservative and the other one on 
conservative flow field variables, are developed for quasi one-dimensional Euler 
equations. In addition to the difference in the representation of the flow field 
variables, the first approach incorporates time integration while the second 
approach neglects the role of time and only takes the converged residual part of 
the solution. These approaches, which are based on the variational methods 
(VM), are used to derive the costate(adjoint) partial differential equations and 
their transversality (boundary) conditions from the differential equations of the 
fluid flow. The costate equations coupled with the flow field equations are solved 
iteratively to get the functional derivative coefficients. Then, these derivative 
coefficients, combined with the flow field variables, are used to find the boundary 
shape which minimizes the performance index (objective functional). 

To demonstrate the method through examples, the shape of the nozzle is 
optimized for the maximum thrust. For this maximization problem, different inlet 
and outlet flow conditions are considered. In the supersonic flow case, the gain in 
thrust is remarkably high. Even in the shock and the subsonic flows, the 
improvement of the thrust is found to be substantial. As demonstrated through 
the cases investigated, a new improvement is that the present variational shape 
optimization approach is capable of resolving flows with shocks. 



71 


6.2 Model Problem 

To demonstrate the versatility of the proposed approach, the time-dependent, 
compressible, quasi one-dimensional Euler equations are chosen as the state 
equations. The corresponding adjoint equations (costate) with their transversality 
(boundary) conditions are derived by variational methods. In this design 
optimization problem, the objective functional is the thrust which is given in Eq. 
(6.7), where the cross-sectional area S is the design variable and the governing 
equations of the model fluid problem are given in Sec. 2.5 which are 

^) + ^) = 0 ( 6 . 1 ) 

dt dx 


djpuS) t <?(p» 2 +P)s_ p dS_^Q 
dt dx dx 


and 


d(P 5g .) t d[{P + pe,)uS] _ Q 
dt dx 


The eqs. (6.1), (6. 2), and (6.3) are mass, momentum and energy, respectively, 
and p, u, e„ and 5, are density, velocity, total energy and cross-sectional area. 

The example problem here is to find S(x) where 0 < x < L to maximize 

J(Q,S) = )){PS}dxdt 
00 


(6.4) 
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subject to 

S(0)<S 0 and S(L)<S, (6.5) 

where L is the dimension of the geometry and T is the maximum time of 
integration. 

6.3 Approach 1 : Perturbed State Equations and Performance Indices in 

NonConservative Variables 

In this approach, we make use of Eqs. (6.1) - (6.3) and the cross-sectional 
area S, which is expressed as 

S(x,x D ) = a-f 6[tanh(0.8x-4)] (6.6) 

A change of the design function S perturbs the flow variables, which in turn gives 
rise to the variation of the functional. Then, the functional derivatives (sensitivity 
coefficients) are obtained by perturbation techniques, for which the primitive 
variables are redefined as 


p = p + 8p u = u+8u e, = e t + 8e, (67) 

and the variation of the boundary is similarly approximated as 

S = S + 8S (6.8) 

where p , u, e t , and S are the nominal values and 8p, 8u , 8e,, and 8S are the 
perturbed values. First, Eqs. (6.7) and (6.8) are inserted into Eqs. (6.1) - (6.4). In 
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the resulting equation, only the first-order terms in the perturbed quantities are 
retained. Recalling that the nominal flow values satisfy the flow equations, the 
perturbed flow equations are obtained as (the overbar signs are dropped for 
convenience) 

d{S8p + p8S ) _ d{uS8p + pS8u + puSS) _ 

^ = r i ( 69 ) 

dt dx 


djuSSp + pSSu + puSS) _ <?[s(u 2 <Sp + 2 pu8u + 8P) + (pu 2 + P)&s] 
dt dx 


dS _ JdS 

+— 8p + p8 — 

dx \dx 


and 

d{Se,8p + pS8e t +pe t 8S) d[uS(e t Sp + p8e, + 8P) + {P + pe,)(S8u + m55)] _ 

Jt = Jx =r 5 

( 6 . 11 ) 


Similarly, the functional J in Eq. (6.7) is perturbed. Because the nominal 
functional values along with the nominal flow variables satisfy the flow equations, 
the perturbed functional is reduced to the following forms for the maximization 
problem: 

8J = ]]{S8P + P8S}dxdt 
00 


( 6 . 12 ) 
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Note that the variation of x is not considered in the previous derivations because 
the length of the nozzle is not changed in the design process. For this approach, 
Eqs. (6.1 ) - (6.3) are put in their vector form as 

M = r - ,6,3) 

dt 


where 

Q = [P’PU’Pe, r 


(6.14) 


r denotes the vector of the spatial derivatives as given in Eqs. (6.9) - (6.1 1), and 
Q is the vector of conserved variables. Then Eq. (6.1 3) can be used to define the 
residual for the nonconservative approach as 


- d(SQ) 

G(x,t,S,Q) = -±—t-r = 0 

dt 


(6.15) 


To eliminate 5P, from Eq. (6.15), one may augment J with the flow equation 
as 

T L 

JAQXS) = J(Q,S) + J J P (x,t)G(x,t,S,Q)dxdt ( 6 . 1 6 ) 

0 0 


where G is given in Eq. (6.15). Then, the variation of Eq. (6.16) is given as 

T L 

SJ a (Q,X,S) = 5J(Q,S ) + Sj J X T (x,t)G(x,t,S,Q)dxdt 

0 0 


(6.17) 
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Because the length of the nozzle is unchanged and G 8X is equal to zero, then 
Eq. (6.17) is simplified to 


T L 

8J a (Q,X,S) = 8J(Q,S) + J J X 7 (x,t)8G(x,t,S,Q)dxdt ( 6 . 1 8 ) 

0 0 


The substitution of the variational expressions from Eqs. (6.9) - (6.12) into Eq. 
(6.18), gives the following equation for the maximization problem. 


T L 

8J a (Q,l’S) = lj\(SSP + P5S)-X l 
00 


d(uS8p + pS8u + pu8S) d(S8p + p8S) 
dx dt 


T L 

'll 


0 0 
T L 


d^Su 2 8p + 2puS8u + SdP + (P+p« 2 )<5sJ d(uS8p + Sp8u + pu8S) 
dx dt 


xdt 


\dxdt 


oo L ^ 


SP + PS ^dx) ^ 


T ( L f \^ [ d[uSe,8p + S(P + pe,)8u + u(P + pe,)8S] d(Se t 8p + pS8e, +pe,8S) 

- iim * + * 


\dxdt 


00 


(6.19) 


The variation of P can be written in terms of the variation of the other flow 
variables with the help of an equation of state, such as the perfect gas law, 

P = p(y-l)[e, -0.5 m 2 ] (6.20) 

Now, Eq. (6.19) can be integrated by parts. For the maximization problem, the 
costate equations, the corresponding boundary (transversality) and initial 
conditions, and their sensitivity coefficients, respectively, are expressed as 



76 


The costate equations are given as 


3 1 + u + e, + u + u 2 ^2. + [u( y^ - o. 5( y - l)u 2 )] 
dr dr dr ax dx 1 J 


dA 


dx 


3 _ 


-(y-l)(e. -0.5u 2 ) 
( 6 . 21 ) 


(l + U )^3. + p^. + 2pu^2- + [(y-l)p i , J -p f ,-P]^i = (r-l)px (6.22) 

dr dxdx dx 


dA 3 dA 3 / v 

dr dx 


(6.23) 


The boundary conditions of Eqs. (6.21) - (6.23) are given as 


r L 

J {sjuA, +u 2 X 2 +(ye, -0.5(y- l)« 3 )A 3 ]} o dr = 0 


(6.24) 


T L 

J {[pSA L + 2puSX 2 +S[ype, -l.5(y-\)pu 2 ]X 3 ^dt = 0 


(6.25) 


J [SX^dt = 0 


(6.26) 


The terminal conditions at time T of Eqs. (6.21 ) - (6.23) are given as 


L 

J + uS&2 + A 3 ]j dx = 0 


(6.27) 


Lf 

J {SX 2 } T dx = 0 


(6.28) 
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L 


| {5A 3 } T <£t = 0 


(6.29) 


The functional derivative coefficients are finally expressed as 


SI ff[ <^1 <?A 2 <?A, (d 2 '\ *^-2 1 52CJ ^ 

•"lir*" x +pe -ir +p “ar- +(p+p " Jarr** 


(6.30a) 


<95 

Because 5 = 5(x,x D ), one obtains 55 = — — 5x 0 , and Eq. (6.30) can be simplified 

ox D 

as 


dJ a r r [ 5A, 5A 2 5A 3 5A, / 2 \ 5A 2 1 dS 

— - = \w^ + P^^ + Pe l -^ + pu-^- L + [P + Pu 2 )— l >- — dxdt 
dx D J 0 J 0 1 5r P dt H ' dt H dx v H ’ dx\dx D 


0 0 
T L 




(6.30b) 


Eqs. (6.21) - (6.29) are solved for an already known flow solution from the 
flow analysis and given boundary and terminal conditions, i.e., Eqs. (6.24) - 
(6.29). From Eq. (6.30), note that the computation of the functional sensitivity 
requires the prior solutions of the state and adjoint (costate) equations. 


6.4 Approach 2: Perturbed State Equations and Performance Indices in 

Conservative Variables 


Unlike the first approach, this approach deals directly with the variation of the 
conservative fluxes with respect to the conservative field variables, lo realize this 
objective, let us recast Eqs. (6.1), (6.2), and (6.3) in a form amenable to 
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conservative formulation. If one follows this procedure, the time dependent quasi 
one-dimensional Euler equations in conservative form are given as 


foe ) t be 

dt dx 


-H= 0 


(6.31) 


where 

E = S[pu, p + pu 2 , (p + pe, )u] T 


(6.32) 


H, = H s (S,Q) = 


f, dS 
0,p—,0 

dx J 


(6.33) 


are the flux and the source terms, respectively. For the derivation of the 

a(se) 

sensitivity, we consider only the steady-state condition, i.e., - = 0. Now, let 

at 

us define a functional, J , as the integrated force along the given contour 5 as 

J = ](PS)dx (6.34) 

0 


Also the residual R is given as 




(6.35) 


Then, with Eqs. (6.34) and (6.35), the functional is modified by 


Lj u 

J a =\ ( PS)dx + J X T Rdx 


(6.36) 


Consequently, the variation of the functional is expressed as 
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5J a = jjps + S(PS) + [(A r + 8Z T ){E X + 5E X -H s - 5H S )] -PS- Z t [E x - 
0 1 

(6.37) 

dE 

where E = — . By neglecting the second-order terms and only keeping the 
1 dx 

linear parts of Eq. (6.37), one recovers the Euler equations in term of the 
variation of the Lagrange multipliers as 


8Z t (E - H ) = 0 in 0 <x<L 


(6.38) 


Then the variation of the objective functional becomes 


SJ a = J |5(P5) + [(Z T ){SE X - 5tf,)]]dx 


(6.39) 


Let us define the following terms for later use: 


dE_ „ dH , 
dQ 


A = ^ B h =- 

h dQ 

dS .. dH s 

s = — M, = — — 

dx ds„ 


5S = ^-8x d 8 

dx 


dS 

dx 


C-£. 

dQ 


9 (dS )K 


dx D \dx J 


(6.40) 


The insertion of Eq. (6.40) into Eq. (6.39) and further simplification of the 
resulting equation leads to 


&I. = 1\S^SQ + PSS + X' 


( AS Q) r ^SQ.j^ SM 


\dx (6.41) 
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Integration by parts of Eq. (6.38) gives 



-X[A - A 7 


BQ 


SQdx+ A 7 /i<52] o +J 


S^8Q + P8S- A 7 ; 


dP 
BQ 


B(*,) 




(6.42) 


where A, B h , C, and M k are matrices given in Appendix C. If one uses Eqs. 
(6.40) and considers only the first variation of the functional in Eq. (6.42) to be 
zero, then one gets the adjoint (costate), boundary conditions and the functional 
sensitivity equations, respectively, as 


-A T X x -B h T i + SC T = 6 (6.43) 

[A 7 A]^=6 (6.44) 


L 

fi/ fl =J[/ > 55-A 7 M,& I ]d[x 

o 


(6.45) 


The insertion of the various partial derivatives of S from Eq. (6.6) into Eq. (6.45) 
gives 


BJ_a_ 

dx D 


=1 


P-^--i T M L ^~ dx 


0L 


dXr 


9Xr 


(6.46) 
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6.5 Solution Algorithm and Surface Modification 

The state and costate equations are discretized using two vector-splitting 
methods: Beam-Warming flux-vector-splitting method for the first approach and 
the Van Leer vector-splitting method for the second approach. In both 
approaches we have used the speed of sound to serve as a sensor to switch 
between different flow regimes at every grid point. 

After adding an appropriate time term X, into Eq. (6.43), the numerical 

integrations of the adjoint equations, i.e Eqs. (6.21) - (6.29) and (6.43) - (6.45), 
are accomplished from the maximum value of L to 0 and from the maximum time 
value T to 0. To bring this in line with the flow equations, the following spatial and 
temporal transformations are carried: 

C = L-X t = T — t (6.47) 


Then, the discretized forms of the adjoint equations are first transformed by the 
above relationships, i.e. Eq. (6.47), and then integrated in time until the 
convergence criteria are met. 

Based upon the steepest descent method, the optimal cross-sectional area 
can be obtained by the recursive relation, which is given as 



where = { Eqn. (6.30b) or (6.46) } 
dx D 


(6.48) 


and e is the step size. 
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6.6 Numerical Results and Discussion 

As the benchmark of the present variational approach, the thrust was 
maximized for three different flow regimes. The first regime is a purely supersonic 
flow with the inflow conditions at M = 1.5. The inflow conditions for the second 
flow regime are identical to those of the first case, but the outflow plane is 
prescribed as a subsonic flow, approximately M = 0.43, to form a shock in the 
duct. The third flow regime is designed to be a purely subsonic flow with inflow 
conditions at M = 0.3. It is observed that the stability limit for both the flow and 
costate analyses depends on the type of inflow and outflow conditions. Secondly, 
even though the flow and the adjoint equations were both hyperbolic partial 
differential equations, their stability criteria were completely different (refer to 
Chap. 5). 

The convergence history for the flow and costate equations are also different. 
In the three cases presented, the rate of convergence depends on the type of 
inflow and outflow conditions imposed at the inlet and exit. Although the time 
integration of the costate equations is dependent on the fluid flow, its quick 
convergence is assured once the flow solution has converged to steady state. 

The thrust evolution for the three cases considered are presented in Figs. 
(6.1) - (6.3). Despite the presence of a strong shock at the middle of the flow 
field, the evolution of the area (Fig. 6.4) and the distributions of the Mach number 
(Fig. 6.5) and pressure (Fig. 6.6) are all smooth. Also in Figs. (6.7), (6.8), and 
(6.9) are presented the area, Mach number, and pressure histories, respectively, 
for the subsonic flow condition. The thrust first increases and then seems to 
decrease and settle down as the optimization progresses. The initial increase 
and then decrease of the thrust for the subsonic case may be attributed to the 
fact that the physics of flow was not accurately modeled for the type of nozzle 
configuration considered. A much better thrust augmentation is achieved for the 
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supersonic flow condition (Fig. 6.1) as compared with the other flow conditions. 
This may be attributed to the noticeable shape change at the upstream location 
in comparison with the downstream location (Fig. 6.10). One also observes that 
the downstream Mach number (Fig. 6.11) increases and the static pressure 
decreases (Fig. 6.12), which results in the increase of the thrust (Fig. 6.1). 

To build further confidence in the present variational shape optimization 
approach, the mass conservation at every cross section of the initial and 
optimized flows were compared (Fig. 6.13). Except for the shock region, the 
mass flow rates were found to be constant. The optimized shape produced a 
sharper and narrower mass jump at the shock area as compared with the initial 
shape. This may suggest that the shock strength attenuation by the optimized 
shape is improved as compared with the initial shape. 

During the numerical experimentation of the costate equations, the boundary 
conditions of the first approach were easier to numerically implement, but the 
boundary conditions of the second approach are superior, especially when there 
is a flow discontinuity (shock) in the flow field. Apart from this apparent 
difference, both approaches give identical results for the optimal solution. 

6.7 Conclusions 

A proof-of-concept study for a new design optimization method based on the 
variational methods has been conducted. The method has been demonstrated 
for the quasi one-dimensional Euler equations. The general design optimization 
incorporates the optimization, CFD analysis, and the adjoint equations analysis. 
The optimization is based on the steepest descent method, and it is intrinsically 
coupled with the flow and adjoint solutions. 

The thrust maximization for purely supersonic, purely subsonic, and mixed 
supersonic-subsonic cases presented demonstrate that the optimized shapes 
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and flow variables are efficiently predicted even with the presence of a strong 
shock in the flow field. They also suggest that the whole optimization needs 
relatively small incremental computer time and memory in addition to the CFD 
analysis and therefore, suggest that the present variation methods is an 
efficient alternative to perform fluid dynamic design optimization for all types of 
flow regimes. 
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Fig. 6.6 Initial and optimized pressure distributions for the shock flow. 
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Fig. 6.7 Area evolution for subsonic flow. 


Pressure 
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Fig. 6.8 Initial and optimized pressure distributions for subsonic flow. 
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Fig. 6.1 1 Initial and optimized Mach number distributions for supersonic flow. 
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Fig. 6.12 Initial and optimized pressure distributions for supersonic flow. 
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Fig. 6.13 Mass error for the shock flow. 
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Chapter 7 

RESULTS AND DISCUSSION ON DESIGN OPTIMIZATION OF INTERNAL 
FLOWS USING TWO-DIMENSIONAL EULER EQUATIONS 

The main thrust of this chapter is to briefly discuss the the numerical results of 
the variational sensitivity analysis that are obtained by the use of two-dimensional 
Euler flow equations. Additionally, the efficiency and accuracy of the variational 
sensitivity in comparison to the finite difference are analyzed. 

7.1 Two-Dimensional Nozzle Optimization Problem Formulation 

At least a couple of reasons can be given for choosing the two-dimensional 
nozzle geometry in order to demonstrate our point of optimization methodology. 
The first one is that one can easily obtain various types of nozzle geometries by 
simply using already known analytical expressions for different flow conditions. 
The second important reason is also the need to develop a scramjet nozzle 
afterbody for the High-Speed Civil Transport. The third one is the need to 
develop efficient wind tunnels with optimal shapes for various experimental wind 
tunnel applications. The optimization problem demonstrated here seeks the 
optimal shape for the maximum thrust in conjunction with the nonreverse flow 
condition at the exit. Hence, the example problem is formulated as the 
maximization of the functional defined by 

j T = jrdT 

r 


(7.1) 
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with the constraint that the static pressure P at the exit assumes a certain 
percentage of the ambient pressure for maximum expansion at that exit lip of 
the solid boundary. Therefore, the constraint is mathematically posed as 


Gj(Q,n) = j 




l-H 


kr<o 


(7.2) 


7.2 Two-Dimensional Nozzle Flow 

The initial geometry for this internal flow configuration is given in Fig. 7.1. It is 
a supersonic nozzle where only half of the physical domain is considered with 
137 x 69 grid points. It is a convex type of geometry with the smallest area at the 
inlet and a diverging afterbody for supersonic expansion. The only aerodynamic 
inequality constraint considered is the criteria on the static pressure at the exit lip 
of the nozzle to reach a certain percentage of the ambient pressure as a 
necessary condition to avoid any reverse flow from underexpansion as the shape 
evolves during the optimization cycle. 

To assess the variational methods for sensitivity analysis, computational 
efficiency and accuracy calculated by variational methods and finite difference 
are compared. One of the obvious limitations with the finite difference is the 
uncertainty to a priori determine the step size that will give reliable sensitivity 
derivatives. The magnitude of the stepsize is dependent on how accurate one 
needs the derivatives to be. If, for instance, one only needs a 10% deviation from 
the assumed exact derivative, then the step size must be under a 10% range of 
the derivative. In our case of computing the sensitivity derivatives using the finite 
difference, we have assigned the step size to be 0.0001 . 

The x component of the design variables (Bezier control points) are a priori 
computed as being spatially invariable, and the variation of the design variables 
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is allowed only in the y direction. This apparent limitation of the design variables 
must not be a hindrance since addition or deletion of any desired design 
variables in the design domain will produce the same result. To verify this claim, 
two sets of design variables, in addition to the assumed optimal number of design 
variables (in this case the optimum is eight design variables), were investigated. 
The first set was performed by increasing the number of design variables by four 
and the second one by decreasing it by four from the optimal number of design 
variables. Here, the optimal number defined as that number of design variables 
which reproduces the closest shape to that of the initial geometry. 

As presented in Table 7.1, the CPU time and memory requirements of 
complete cycle of optimization for the two additional sets of design variables are 
almost identical for the two-dimensional optimization case. Therefore, the eight 
design variables are considered as the optimal number of design variables which 
produced the desired computational efficiency for our test case. On the other 
hand, this slight memory increase as the number of design variables increases 
could be a warning to the eventual computational memory increase as the 
dimensionality, number of constraints, and design variables increase. The 
second aspect of the role played by the number of design variables may be the 
influence on the optimal shape (Fig. 7.2). All three categories of the design 
variables produced slightly different optimum shapes from each other. Comparing 
all three shapes (Fig. 7.2), the shape produced by twelve design variables 
appears to follow the shape produced by the four design variables in the 
compression area (upstream) and the shape of the eight design variables in the 
expansion area (downstream). The shape generated by the optimal number of 
eight design variables shows a slight change of shape upstream, from 
approximately x=0.1 to x=0.375, and downstream, from approximately x = 0.7 to 
x = 1.0 as compared with the shapes generated by the other sets of design 
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variables. The shape change in the compression area seems to be more 
desirable beause it produces high-pressure ratios and thereby gives more thrust 
as one integrates the change of pressure along the changing nozzle shape. The 
shape change in the expansion region, on the other hand, reduces the ratio of 
the static pressure to the ambient pressure, which results in less thrust 
augmentation. This physical phenomena is further reflected in Fig. 7.3 where the 
optimal thrust of the eighth design variable shape is higher than the other two 
design variable shapes. 

From the parametric studies (four, eight, and twelve design variables) 
conducted, one may conclude that the eight design variables are the optimal 
number of design variables to sufficiently represent the nozzle shape and at the 
same time to give a better thrust and computational efficiency. 

The evolution of the design variables for the variational methods and the finite 
difference approach are given in Figs. 7.4 and 7.5, respectively. Except for the 
second and the seventh design variables, the general trend of the evolution of 
the design variables in both approaches is similar. In the variational case, the 
second design variable approaches the first design variable and the seventh one 
tends to come close to the eighth design variable. In the finite difference case, 
however, the second and the seventh design variables tend to pull away from the 
first and eighth design variables, respectively. As shown in Fig. 7.6, due to the 
movement of the second and the seventh design variables in the opposite 
direction, the optimal shapes of the variational methods and finite difference are 
slightly different. As explained in the parametric studies, the decrease of the 
optimal (as compared with the initial) shape or optimal design variables in the 
compression region is much more advantageous to the decrease of the optimal 
shape or optimal design variables in the expansion region for the supersonic flow 
regime. This is due to the effect that the decrease of the shape in the upstream 
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results in the substantial gain of high pressure ratio (compare Figs. 7.7 and 7.8) 
which favors the augmentation of more thrust (Fig. 7.9) in the design process. 
Figure 7.9 also clearly indicates that the pressure distribution in the expansion 
region in general and at the lip of the nozzle in particular is within the constraint 
specification as imposed in the aerodynamical constraint given by Eq. (7.2). 

As given in Table 7.2, the accuracy of the variational methods is verified by 
comparing the variational functional sensitivity derivatives to the functional 
derivatives of finite difference. If one takes into consideration that the sensitivity 
coefficients of the finite difference are dependent on the step sizes, then the 
gradient values obtained by the variational methods are well within the 
engineering prediction range, except for the second and the seventh sensitivity 
coefficients. The discrepancy of those two sensitivity values may be associated 
with the difficulty to properly implement the boundary conditions of the adjoint 
equations. Despite the differences on these two sensitivity derivatives which 
correspond to the second and seventh design variables, the optimal shape and 
thrust of the variational methods are comparable with those of the finite 
difference as presented in Table 7.3 and Fig. 7.6. It is known that the finite 
difference uses function evaluations to compute the gradient information while 
the variational methods solve another set of partial differential equations and 
sensitivity derivative equations. Due to this, there is a memory increment of 
approximately 1.3 mega words as shown in Table 7.4. This slight increment in 
memory is negligible as compared with the other gradient-based sensitivity 
analysis approaches, such as the discrete sensitivity analysis which requires 
higher memory allocation for the given optimization problem. 



103 


Table 7.1 . CPU Time and Memory for Four, Eight, and Twelve Design Variables 
With Variational Methods 


Design variables 

CPU time (sec) 


4 

868.0463 

5.249459 

8 

864.2226 

5.249939 

12 

866.2128 

5.250579 


Table 7.2. Sensitivity Derivatives by Variational Methods and Finite 
Difference 



Variational 

methods 

Finite 

difference 

Deviation (%) 

1 

9.1483E-2 

9.4441 E-2 

3.1 

2 

7.9228E-2 

1.1062E-2 

86.0 

3 

-6.6563E-2 

-4.7906E-2 

28.1 

4 

-5.5491 E-2 

-5.9409E-2 

6.6 

5 

-4.6421 E-2 

-5.3278E-2 

12.9 

6 

-3.8979E-2 

-4.6287E-2 

15.8 






-3.21 86E-2 


-6.9988E-2 


53.9 
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Table 7.3. Initial and Optimal Values of Functional and Constraint for 
Variational Methods and Finite Difference. 



Variational 

methods 

Finite 

difference 

Initial 

Functional 

0.045481 

0.045481 


Constraint 

- 2.10787 

-2.10787 

Optimum 

Functional 

0.049958 

0.49885 


Constraint 

-0.5858 

-0.5668 


Table 7.4. Efficiency Comparison Between Variational Methods and Finite 
Difference. 



Variational 

methods 

Finite 

difference 


Complete 


865.098 

4356.33 


optimization 




CPU 

Single 

Euler 


128.23 

time (sec) 

analysis 

equations 





Co-state 

58.59 




equations 



Memory 

Complete 


5.25 (with 

3.98 (no 

(MGW) 

optimization 


sensitivity 

sensitivity 




eqs.) 

e q s ) 
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Fig. 7.3 Effect of the number of design variables on the optimal thrust. 
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Nozzle length 


Fig.7.6 Optimal shapes by variational methods and finite difference. 
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Fig. 7.7 Optimal pressure countour by variational methods. 
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Fig. 7.8 Pressure distributions on the solid and symmetry line by variational methods. 
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Fig. 7.9 Optimal thrust history by variational methods and finite difference. 
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Chapter 8 

CONCLUSIONS AND RECOMMENDATIONS 


A versatile design optimization approach must be independent of: (1) the level 
of approximations, methods of numerical integration, and discretization of the 
flow analysis, (2) the grid points and grid sensitivities, (3) the initial design points 
and solutions, and (4) the flow regimes. In the case considered, starting from the 
converged solution of the two-dimensional Euler flow equations, the sensitivity 
derivatives for the optimal solution are derived (See Chap. 3.) The distinctive 
feature of the variational approach is that the converged solution can be obtained 
from any level of approximations, methods of numerical integration, and variety of 
discretizations of the fluid flow. For instance, either the two-dimensional full 
potential or the two-dimensional Euler equations can be chosen for the level of 
the flow field approximation. One can also choose either the Beam-Warming, or 
Van Leer flux-splitting, or Roe flux-differencing, and any implicit or the explicit 
Runge-Kutta time-advancing method can be adopted. The only requirement in 
this approach is that the state solution has to completely converge so that the 
Jacobian matrices that are needed in the costate, transversality, and sensitivity 
derivative equations are accurately constructed. Even though the Jacobian 
matrices in the costate equations are transposed, they are the by-product of the 
converged solution from the state equations, and there is no extraneous effort 
involved in obtaining them. 

The time integration of the costate equations, as described in Chap. 4, 
requires an understanding of the stability behavior of the costate equations. 
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Because of the negative entries of the transposed Jacobians, the Von Neumann 
stability analysis indicates that the costate equations are unstable for the CFL 
and wave number considered. To make the costate equations stable, the 
transformation of the costate's computational domain in line with the state's 
computational domain was undertaken. This unique approach produced a stable 
costate solution that is necessary in obtaining the sensitivity derivatives of the 
functional and constraints with respect to the Bezier control points. 

Variational methods in sensitivity analysis and shape optimization for 
aerodynamic applications have been considered. The main components of the 
approach are the control theory and calculus of variations where the design 
control is the continuous domain boundary represented by the Bezier control 
points or design variables, Eqs. (3.5) - (3.7). Because the domain boundary is 
continuously evolving until it reaches its optimal shape, the contribution of the 
changing domain and domain boundary are used and incorporated in the 
derivation of the optimality conditions and sensitivity derivatives through the 
surface transformation matrix J s (Appendix A) and curvatures (Eqs. (3.24 and 

3.25)). The use of the steady-state solutions of the two-dimensional Euler 
equations results in the derivation of the costate or the adjoint equations, the 
boundary or transversality conditions, and the sensitivity derivatives for the 
generic design optimization problem where the functional and the constraint are 
defined on the domain boundary. The final numerical computation of the 
sensitivity derivatives is carried out for the converged state and costate 
equations. 

Since the stable and converged solution of the costate equations is imperative 
for the successful computation of the sensitivity derivatives and design 
optimization, a complete stability analysis of the Euler and costate equations is 
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treated. Based on the findings of the stability criteria, the stable and converged 
numerical integration of the costate equations is guaranteed only if the domain of 
integration is transformed in line with the state equations or any other numerical 
technique which handles the reverse nature of costate equations is adopted. 

A two-dimensional nozzle optimization problem was considered, and the 
application of variational methods to compute the optimal shape for the maximum 
thrust is presented. During the design process, the supersonic nozzle remained 
supersonic while improving the performance index or thrust (Table 7.2). Also, 
while the VM's computational accuracy (Table 7.3) is comparable with the finite 
difference, its computational efficiency and memory savings (Table 7.4) are found 
to be substantial. As memory and computational efficiency are the bottle-necks 
for large two-dimensional and three-dimensional problem in general, variational 
methods are one of the most viable candidates in solving design optimization 
problems. 

Design optimization requires gradient information of both the functional and 
constraints. But as the number of constraints increases, the computational 
intensity and memory may be prohibitive to do any realistic optimization. This can 
be overcome by converting the constrained problem into an unconstrained 
problem through the introduction of penality function methods. This is highly 
important especially for multidisciplinary optimiztion where the number of 
constraints and design variables is high. Therefore, procedures similar to the 
penalty function methods must be explored to efficiently use variational methods 
in obtaining sensitivity derivatves. To further build confidence in the proposed 
approach, two-dimensional internal and external flows with shocks must be 
investigated. Finally, the concept must be extended to two-dimensional viscous 
and three-dimensional Euler flows. 
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APPENDIX A 

CONTINUOUS DOMAIN IN VARIATIONAL METHODS 

Consider the domain Q and the solution vectors Q which are transformed 
from an original domain Q and from the nominal solution vector Q by the 
following one-to-one mappings [84]: 


X =o(z ; e)=i + e 


d<S> 

de 




c=0 


(A. 1 ) 


X = X + e 


d<& 

de 


0(e') 


£-0 


(A-2) 


where e is a small quantity. Likewise, 


Q = W;e) = TO0) + e 


de 


0(e 2 ) 


e=0 


(A.3) 


Q =Q + e 


dW 

de 


0(e>) 


Je=0 


(A. 4) 


Now, take only the linear parts of Eqs. (A.2) and (A.4), and one gets 


AX = X - X = e\ 


d<& 

de 


= SX 


(A.5) 


Jc=0 
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AG=e(x)-e(x)=/^l «5g 

L j£ = 0 


(A.6) 


The variation of the solution vector in the usual variational sense (i.e., 
variation at the same coordinate location) can also be defined as 


AQ=Q(X)-Q(X) = 8Q 


(A.7) 


Then for the relation between SQ and SQ, the following expression can be 
established as 


aq=S(z)-q(x) 

= d(^]S(X) + U(X)-Q(X) 

=(Q(z)-U(X)) + (§(X)-Q(X)) 


(A. 8) 


From the Eq. (A.8), q[x^-Q(X) can be approximated by 


e(x)-e(*) -§<*-*)-§* 


(A. 9) 


and with Eq. (A.9), Eq. (A.8) becomes 


AQ = ^8X + 8Q = SQ 
aX 


(A. 10) 


Then with Eq. (A. 10), SQ is approximated as 


SQ=SQ-^SX 

oX 


(A. 11) 
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Equation (A. 10) explains that the total variation of the solution vector 5Q is 

composed of parts due to the variation of the state vector 8Q and the coordinate 
dO - 

transformation —8X. 

dX 



131 


APPENDIX B 

DERIVATION OF THE SURFACE TRANSFORMATION JACOBIAN 


For a small quantity e, the transformed space is represented by Eq.(A.I). 
Then by the chain rule, the surface Jacobian J s can be obtained as [84] 


dX ,. . 

— - = 5 + £ 

f 3 [ 

’(d*Y 

A 

dx j 

l dX j ■ 


e-o J 


J = 1.2, .... m; / = 1, 2, .... n (B.1) 


where S t j is the Kronecker delta equal to 1 if i = j and 0 otherwise. The second 

term on the right hand-side of Eq. (B.1) indicates the contribution due to the 
space transformation. Then, the total change of the area due to the space 
transformation is the determinant of Eq. (B.1 ) which is expressed as 






f d [ 




1 

l®il 

UeJ 

j _ 

£=0 J 

1 


j =1,2, .... m 


(B.2) 


Further simplification of Eq. (B.2) gives 


f + ef< 

s ( 



)] 

i=] 

1«A 

.1 de j 

i _ 

e=oJJ 


(B.3) 


where I is the identity diagonal matrix, because the variation of the boundary is 
due to the variations of the boundary coordinates and that 



then, Eq. (B.3) can be expressed as 


\Js 


I + V.5X 


(B.5) 


which is the total change of the surface due to the change of the variable domain. 
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APPENDIX C 

VARIOUS MATRICES IN CONSERVATIVE FORMS FROM THE QUASI ONE- 
DIMENSIONAL EULER FORMULATION. 


The various matrices defined for the quasi one-dimensional Euler equations in 
the Cartesian coordinate systems are asfollows: 

1. By use of Beam-Warming flux splitting, the approximate Jacobian matrix A is 
computed as [76] 


A = S\ 


0 


1 

-(y-3)u 


, 3 3(7-1) 2 

(y-\)u 3 -yue, ye, — u z 


0 

(7-1) | 

yu 


(C.1) 


2. By use of the Van Leer flux splitting, the exact Jacobian matrix A is computed 
as [80] 


a = local speed of sound 

E = [^1,VL>^2,VL’^3,Vl] 

is the flux vector where the components are defined as 


p{u + pm-af _ {q 2 + pm*aq l f 

^\ % vl ~ pm ' A ~ pm * A 

4 a 4 aq } 


4 a 


(C.2) 
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^2 ,VX — E\ % VL * 


[(y-l)u + 2 pm a] „ [{y-l)q 2 + 2pmaq l ] /r% ^ 

— t l VL ■ 

r w, 


e- [(r-l)» + 2'^-flf_ r [{y-l)q 2 + 2-pm aq.f /f> ^ 

3VL UVL ' 2{f-\) ~ hVL ' 2(y 2 -\)qf ( } 


where the parameter to switch between positive and negative contributions is 


pm = ( ± ) 

{0.5) 

and Q and P are defined as 


Q = [p,pu,pe,f =[<7i,<? 2 ,<73 f 

(C.6) 

and 

P = (Y-l)\pe,-p~] = (Y-\) 

L 2 J 2c h\ 

(C.7) 


The derivatives of the fluxes with respect to the conservative variables are 


3E\,VL _ pm 

~da~~ 


' fa+pm-aqjdiqt + pm-aq,) _ ( )2 

aq { dQ r W2 da 


(C.8) 


d^i,vL _ (/ 2 • ptn • dEyy^ 


dQr 


M 


dQr 


+ 


E U vL\d({y-l)q 2 + 2 -Pm ^) „ \dq, 

— — --((r-i)<f 2 + 2 P ma( 7i)- 


Q\dQ r 


dQ, 


(C.9) 
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,VL __ ( Y ^ ' pm ' QQx ] ^E\,VL , 

2(y 2 -l) ? , ! 5 gT 


£ i.w. f 2 f(i'~ 1 ) i ?2 +2 P m ' ai ?] £(jQ | 
2(y ! -i)l «. 2 <*2, J 


Then the Jacobian components are determined from 


SQ, 


dE\yL ^^i,VL 

dQr ’ <5Qr ’ dQ r 


(C.10) 


(C.11) 


where r is the index looping from 1 to 3 and the other quantities are defined as 


*„ = (r-iK 




0 

Qi 

Qx 

0 


0 

1 

0 


C = (y- 




(C.12) 


(C.13) 


M* = [0,/>,0f 


(C.14) 


where P is given in Eq. (C.7). 
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APPENDIX D 

TWO-DIMENSIONAL JACOBIAN MATICES THAT ARE USED FOR THE 

STABILITY ANALYSIS 

By use of the Van Leer flux-splitting technique in Cartesian coordinates, the 
exact Jacobian matrices A and B of the inviscid two-dimensional Euler Equations 
are computed. To realize this, let us define the field variables, fluxes, and speed 
of sound in the conservative form as 


q = [p.p^pe.r = 


(D.1) 


a = < 


7(7 ~ 1 ) 


tft 

tfi 


/■ 2 X 2^ 
<72 +<73 

v. 2 <?i , 


Yl Vi 


(local speed of sound) 


(D.2) 


M = — , (Mach number) 
a 


(D.3) 


where the parameter to switch between positive and negative contributions is 
pm = (±) (D.4) 

and 




- J — + pm 

IP * ) ... 


M. 


pmq x a • 


— - + pm 
Ui 


(D.5) 
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_f (y~ 1 K + 2 -pm-a ~l_r (r- l )g; + 2 ■ pm - aq : 


(D.6) 


/£— = fL.fi (momentum components) 


(D.7) 


(f± \ =f ± \ [(r-l)Uj + 2-pm a] 2 | f yM _ 

/*mi 2(y 2 -l) V^J 

[(y-l)(7: + 2'p/ri’fl] ^ 

/—, 1 r7~2 7T 2 — L+ T2 (energy in the x direction) (D.8) 

2(r -1)<7: (,2(7, J 


, s , \(r-nu, + 2pmaf (V'j 

\J'~r t y) y Jnuvs 2(f ~ l) ^ 2 J 


(D.9) 


[(r-i)<7> + 2 p/w a] 2 

2(r 2 -iK + U«i 


(energy in the y direction) (D.10) 


Then the flux vectors in the x and y directions, respectively, are given as 


f = r f ± ^ ft i 7 

VL I •/ rruiss ’ J energy I 


(D.1 1 ) 


The derivatives of the fluxes are 

dfL* P m f Y^[^r fl ] n./ \( Mj ^ , <7, P 

dQ r 4 [q, v J dQ, Wl \ qi p J dQ r 


(D.12) 


fttomj _ f ± tfL* , f ± 


(D.13) 



138 


d(flj. 

[(y-VQj+l- Pmaf 

f q * 1 

dfL, 

dQ r 

2(f- IK + 

wJ 

dQr 


• + 


/; 


± 

mass 


[( 7 - l )? y +2 pm of 

f 2 \ 

2 ( r 2 - i )?! 2 

J 


dQ, 




L - 


dQ, 

d 


[(7“!)?; +2-pm*fl] 

( ^ l 


2(^-1)?, 2 + l 

wJ 



f\ 


± 

mass 


[(7 — 1)^ + 2- pm • tz] 2 

f ?2 

2(y 2 - 1 K 



dQr 


(D.14) 


(D.15) 


Then, the A and B Jacobian matrices in the x and y directions, respectively 
are constructed as 


dFyL,x _ 

df mass },mom Of energy, x 

dQr 

_ dQ, ' dQ, ' fa J 

dFy\, _ 

mass $ J,rnom energy, y 

dQr 

. dQ, ' dQ, ’ <?a J 


(D.16) 


(D.17) 


where r is the index looping from 1 to 4, and j from 2 to 3. 
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APPENDIX E 


DERIVATION OF TWO-DIMENSIONAL CONSERVATIVE VAN LEER 
JACOBINAN MATRICES IN GENERALIZED COORDINATES 


Hyperbolic equations such as the Euler equations of the fluid flow, are 
numerically treated based on the direction of the physical propagation of waves. 
Therefore, the fluxes and the Jacobian are split into the positive and negative 
contributions and are discretized according to the sign of the Jacobian matrices. 
The Van Leer flux-vector-splitting technique in generalized coordinates, which 
follows this upwinding method and which has the smooth differentiability property 
at sonic transitions and stagnation points, is used for our flow analysis and co- 
state equations. For the general derivation of the various terms in the splitting 
procedure, we define the following quantities: 


* |Vfl 
(E.1) 


(direction cosines, k = x , y ) 


rlen = + £j) (cell area) (E.2) 

U = + = (^ x + (contravelocity) (E.3) 

P Qx 


Q = [p,puj,pe,] T = [q„q 2 ,q„q 4 ] T (conservative variables) 


(E.4) 
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a = < 


7(7 - 1 .) 


<h 

<h 




K 


2q 


i yj 


(local speed of sound) 


(E.5) 

M, = — , Mach number (E.6) 

4 a 


where the parameter to switch between positive and negative contributions is 
pm = ( ± ) (E.7) 

(m c + pm] (M c + pm) 

fLs =rlen pm- p a - — — = rlen- pm- q x - a- (E.8) 



—i-U + 2 a ■ pm) + — 
7 


(E.9) 


/* = /* /* 
^ j , mom J massJ ] 


(E.10) 


/* 

v en 


= f ± 

energy J mass 


(1 - y)U 2 + 2 • pm ■ (y - 1 )aU + 2a 2 

( qI + qI \ 


l 2<? 2 . J_ 


(E.1 1) 


\ = 1 V C| f r± r± r± 1 T 

1 J [y mass'* j j , mom* J energy J 


(E.1 2) 


dMj _ 1 

*21 a 2 


dU da 
a- U- 


dQ r dQr 


(E.13) 


By the use of Eqs. (E.1) - (E.13), the derivatives of the mass, momentum, and 
energy with respect to the conservative field variables can be expressed as 
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= rlen * — 
dQr 4 


a ( M { + || + K + H 2 H + 2aq ■ K + H 


\2 <?a 




(E.14) 


dQr 


c» 


rL 32- 


dU da 

+ !■ pm 


dQr 


+ * 






(E.15) 


_ f ± <yLx , f * dg_ 

dQ r jj dQ r Jmaa dQ r 


(F.16) 


2 = M^ll 

2<7| 2 


(E.17) 


df energy _ /-± df mass ± [ (a • JVn U) dU_ ^ 

^ W o>Q r / H (y + i) <?a 


U + 2 


1 <?a dq22 


{y-\)\{y + \)dQ r ' dQ r 
(E.18) 


With Eqs. (E.14) - (E.18), the general expressions for the split Jacobian in 
conservative variables can be symbolically represented as 


dFvL_ 


dQr 


dfLs dj % mom df, 


energy 


dQ r ’ dQ r ' dQ r 


(E.20) 


where r is the index looping from 1 to 4 and j from 2 to 3. When the computation 

dF ± 

is performed in the E, dierction, — denotes the split Jacobian matrices in the 

dQr 
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£ direction, similarly, when the computation is performed in the r\ direction, 

— — denotes the split Jacobians in the ^direction. 

9Q r 



